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1.  INTRODUCTION 

In  recent  years  there  has  been  an  Increasing  interest  In  finite  element 
software  evaluations  due  to  the  widespread  availability  of  such  packages. 

Also  recognizing  this  need,  an  ISEG  interagency  Software  ^valuation  Group) 
associated  with  the  armed  forces  and  other  government  organizations  -  such  as 
the  National  Science  Foundation,  the  Department  of  Energy,  and  the  Nuclear 
Regulatory  Commission  -  was  formed  for  continuous  coordination  of  evaluating 
engineering  application  software.  In  this  connection,  a  set  of  criteria  for 
selecting  software,  selecting  contractors  and  evaluation  procedures  were 
developed  by  Nickel  1  [1,2].  Subsequently,  three  general  purpose  finite 
element  programs  for  structural  mechanics  applications  were  selected  and 
evaluated.  These  are  ADINA,  NASTRAN  and  STA6SC-1  and  the  evaluation  results 
were  reported  in  references  [3,4,5].  A  comprehensive  summary  of  evaluation 
work  is  given  by  Nickel!  [2]. 

As  a  part  of  ISEG's  continuing  effort,  AiiAQUS  was  also  included  for 
evaluation.  ABAQUS  is  a  general  purpose  finite  element  code  for  nonlinear 
static  and  dynamic  analysis  of  structural  problems,  as  well  as  steady-state 
and  transient  analysis  of  heat  transfer  problems.  The  code  was  a  relatively 
new  entry  to  the  nonlinear  structural  software  market;  it  was  developed  and 
released  by  Hibbitt  X  Karlsson,  Inc.,  1979.  Although  the  version  3-  ABAQUS* 
has  limited  analysis  capabilities  (as  compared  to  MARC  [6]  for  example),  it 
has  the  potential  to  become  a  major  nonlinear  code  due  to  its  modern  concept 
software  design  and  strong  commitment  in  coding  development  by  the  company. 


*It  is  noted  that  only  the  version  3  of  ABAQUS  was  available  for  the 
evaluation  work.  Therefore,  whenever  ABAQUS  is  mentioned  again  in  this 
report,  it  implies  the  version  3  of  ABAQUS. 
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By  following  the  procedure  outlined  in  [1],  the  evaluation  work  consists 
of  the  following: 

i)  Review  of  documentation 

ii)  Program  organization  and  data  base  design 

iii)  Functional  description 

iv)  Benchmark  runs  for  coding  verification 

v)  Benchmark  runs  to  determine  the  numerical  characteristics  of  ABAQUS 
(advanced  evaluation). 

Since  ABAQUS  ’?)  marketed  as  a  production  code,  unlike  ADINA  primarily  us  a 
research  code,  no  attempt  was  made  in  the  review  of  detailed  program 
architecture.  Instead,  discussion  is  focused  on  the  overall  coding  structure 
and  its  data  base  design.  This  is  in  contrast  to  our  previous  review  work  on 
AUINA  [3].  Tne  purpose  of  functional  description  is  to  give  a  comprehensive 
overview  of  the  solution  methodologies  adopted  and  analysis  features  available 
in  ABAQUS.  The  section  on  verification  problems  is  to  demonstrate  the 
operational  status  of  double  precision  version  of  ABAQUS  on  the  IBM  main 
frame.  Finally,  a  great  deal  of  effort  was  devoted  to  advanced  evaluation  in 
which  a  selected  number  of  benchmark  problems  were  run  to  determine  the 
operating  characteristics  of  ABAQUS  in  terms  of  its  solution  convergence 
behavior  and  computational  efficiency  (i.e.  CPU  time). 
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2.  Documentation 

ABAQUS  is  supported  by  four  documents: 

•  User's  Manual  [7] 

•  Theoretical  Manual  [8] 

•  System  Manual  [9] 

•  Example  Problem  Manual  [10] 

2.1  User's  Manual 

Obviously,  the  purpose  of  user's  manual  is  to  provide  specific 
instructions  on  data  input  and  to  give  a  description  of  output  features.  The 
major  importance  of  this  manual  is  its  clarity  in  instructions. 

The  user's  manual  contains  ten  sections.  Sections  1-4  are  the  reference 
sections  in  which  input  organization  and  syntax  rules,  analysis  procedures, 
element  and  material  libraries  are  described.  Sections  6-10  contain  data 
input  instructions.  Section  5  is  not  included  in  the  manual;  it  may  have  been 
reserved  to  cover  additional  features  which  will  be  added  into  the  program. 

For  the  most  part,  the  manual  was  clearly  written  and  easy  to  follow.  To 
prepare  the  input  data  for  a  problem,  the  user  should  be  able  to  find 
sufficient  information  from  the  instruction  sections.  To  input  element  and 
material  data,  some  cross-reference  in  Sections  3  and  4  is  often  necessary. 
Although  the  manual  is  generally  well  documented,  two  points  may  deserve  some 
future  attention: 

i)  General  shell  section  (6.5.9)  is  provided  to  allow  direct  input  of  a 
shell  cross-section  with  elastic  anisotropic,  layered  or  corrugated 
material.  However,  the  definition  of  section  stiffness  matrix  D  was  not 
clearly  specified.  Confusion  exists  in  deciding  whether  the  D-matrix 
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for  larje  deformation  analysis  has  the  same  meaning  as  the  one  for  small 
deformation.  An  example  is  needed  in  order  to  clarify  this  point, 
ii)  Restart  (6.8.4)  option  can  be  exercised  to  save  and  reuse  data  and 
analysis  results.  The  manual  does  not  clearly  illustrate  how  to  use 
♦RESTART  option  to  conduct  a  multi-step  analysis.  For  example,  in  the 
case  of  dynamic  analysis  it  is  not  clear  that  whether  the  analysis  time 
should  be  accounted  from  the  ve'y  beginning  of  the  analysis  or  from  the 
starting  point  of  each  new  step?  Similar  question  arises  for  the 
definition  of  new  amplitude. 

2.2  Theoretical  Manual 

The  theoretical  manual  was  written  for  the  advanced  user  who  has  had 
sufficient  background  in  (nonlinear)  continuum  mechanics  and  wishes  to  know 
the  theory  behind  the  code.  The  manual  include  the  following  basic  items; 

•  Review  of  basic  mechanics  principles 

•  Description  of  analysis  procedures 

•  Element  formulations 

•  Constitutive  theories 

•  General  formulations  for  heat  transfer  analysis 

All  basic  assumptions  are  properly  stated.  Kinematic  and  kinetic  variables 
are  clearly  defined.  This  is  a  very  important  point,  especially  for  large 
deformation  analysis. 

In  the  element  formulations,  considerable  details  are  given  for  the  beam, 
shell  and  elbow  elements.  Discussion  on  solid  elements  is  rather  sketchy; 
related  literature  should  be  cited  to  compensate  this  effect. 

2.3  System  Manual 

The  system  manual  contains  a  fairly  detailed  description  on  the  program 
design  of  ABAQUS,  The  purpose  of  this  manual,  as  stated  in  [9],  is  intended 
for  a  user  to  have  access  to  the  source  code.  The  user  may  be  a  system 


2-3 


proyratnmer  who  is  charyed  with  inst.t11:;ij  program  at  a  computer  site,  or  a 
person  who  will  make  changes  or  additions  to  the  code  in  order  to  solve  a 
particular  class  of  problems.  The  manual  contains  the  following  important 
information : 

•  Overall  program  design  concept  and  decisions 

•  Brief  review  of  theoretical  basis 

•  Detailed  flow  of  the  PRE  and  MAIN  Programs 

•  Description  of  data  base  design  and  data  base  managers 

•  Logic  for  allocation  of  workspace 

•  Description  of  labelled  common  blocks  in  PRE  and  MAIN 

•  Definitions  of  subroutines  used 

•  User  interfaces 

While  the  reviewer  found  that  the  system  manual  was  helpful  to  unde' 
certain  part  of  the  source  code,  a  good  portion  of  the  manual  has  become  out¬ 
dated.  A  number  of  subroutines,  which  were  added  to  or  deleted  from  the 
program,  are  not  mentioned.  Therefore,  the  system  manual,  as  the  way  it 
stands,  can  at  most  be  used  as  a  reference  material  until  it  is  properly 
updated. 

2.4  Example  Problem  Manual 

An  example  problem  manual  which  contains  about  sixteen  sample  problems  is 
also  available.  The  purpose  of  this  manual,  as  it  appears,  is  to  assist  the 
new  user  in  setting  up  problems,  or  to  verify  the  code  upon  installation. 

Each  problem  is  presented  in  a  standard  format:  problem  statement; 
description  of  material,  geometry,  loading  and  boundary  conditions;  results 
and  discussion;  listing  of  input.  In  general,  the  example  problems  are 
clearly  explained  in  the  manual.  However,  most  of  the  problems  are 
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concentrated  on  the  use  of  beam  or  pipe  elements  with  elastic  or 
elasti-plastic  materials.  There  is  a  lack  of  representation  of  different 
elements  or  material  models.  Furthermore,  the  manual  that  the  reviewer  has  is 
in  a  rough  draft  form  which  does  not  follow  a  proper  sequence  order. 


3.  PROGRAM  ORGANIZATION  AND  DATA  BASE  DESIGN 


ABAyuS  was  designed  primarily  as  a  production  code  for  solving 
nonlinear  problems.  The  program  must  be  run  on  a  fairly  big  computer  (such  as 
a  mainframe  or  a  super-mini  computer)  so  that  a  reasonably  large  core  size  can 
be  allocated  in  order  to  maintain  its  efficiency.  In  terms  of  its  coding 
structure,  A8AQUS  is  divided  distinctively  into  two  sub-programs:  namely  PRE 
and  MAIN  programs.  The  major  functions  of  PRE  are:  i)  to  read  the  input  data 
and  sort  the  data  onto  a  data  structure  from  which  the  fiAIN  will  read,  and  ii) 
to  perforin  data  check  analysis.  On  the  other  hand,  the  function  of  the  MAIN 
program  is  strictly  for  engineering  analysis.  The  obvious  advantage  of  this 
division  is  that  the  PRE  can  be  repeatedly  run  alone  for  data  checking  with 
quick  computer  turn  around  time.  This  feature  is  extremely  helpful  for 
computer  systems  which  assign  higher  priorities  to  jobs  with  short  CPU  time 
and/or  small  central  core  usage. 

In  addition,  A8AQIJS  utilizes  an  extensive  data  base  design  concept  in 
the  development  of  its  coding.  This  approach  provides  three  apparent 
advantages:  i)  protection  for  data  file  over-write,  ii)  efficient  utilization 
of  fast  core  storage,  and  iii)  ease  of  future  coding  extension  or 
modifications.  In  this  section,  we  will  review  each  of  the  following  three 
items  in  ABAQUS;  i)  PRE  program,  ii)  MAIN  program,  and  iii)  data  base  design. 

3.1  PKE  Program 

The  major  functions  of  the  PRE,  as  stated  in  the  ABAQUS  system  manual 
[9],  are  to  perform  data  check  and  create  the  data  bases  for  the  MAIN.  The 
data  bases  are  written  onto  a  communication  data  base  which  are  then  passed  on 
to  the  MAIN.  If  any  fatal  errors  are  detected,  job  execution  will  be 
terminated  at  the  end  of  PRE. 
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The  PRE  proijrdiii  consists  of  9-moclu1es  which  are  executed  sequentially. 
These  modules  are  differentiated  by  step  designations: 

Step  1.  Scan  input  data  -  The  PRE  reads  all  input  data  and  provides  an 
'echo'  printout  of  user's  input  data.  '11  the  common  blocks  are  initialized 
and  the  locations  of  the  data  bases  are  defined.  Each  data  base  has  a 
constant  length  of  1024  words  except  that  the  node  point  data  base  has  a 
length  of  2048  words  (double  precision  for  IBM  computer).  Also,  the  buffers 
of  the  data  bases  are  initialized. 

Step  2.  Read  bulk  data  -The  input  data  are  read  again  and  they  are  separated 
into  two  different  types:  a)  Bulk  data  -data  associated  with  elements,  nodes, 
element  sets,  and  node  sets,  and  b  )  short  list  data  -data  other  than  a.  Mesh 
and  incremental  generations,  if  encountered,  are  performed.  Information  such 
as  BEAM  SECTION,  MATERIAL,  STEP,  etc.  are  written  onto  an  input  data  base,  and 
an  input  data  manager  is  also  created.  The  buffer  space  for  the  input  data 
base  manager  was  set  to  1500  words  in  the  program;  3  -words  required  for  each 
input  keyword  card.  Therefore,  a  maximum  of  500  keyword  cards  can  be  read  in 
for  the  current  version  of  ABAQUS  (version  3)  unless  the  buffer  space  is 
enlarged.  The  short  list  data  is  separated  into  different  volumes.  At  the 
end  of  this  .tep,  volume  1  of  short  list  data  base  is  extracted  from  a  data 
block  in  subroutine  BELTfP,  and  then  copied  onto  the  buffer  area  of  the  short 
1 ist  data  base. 

Step  3.  Read  material  properties  -Material  property  data  are  read  and  written 
onto  volume  :i  of  the  short  list  data  base.  Volume  3  of  short  list  data  base 
is  also  defined. 

Step  4.  Build  element  data  base  -  Element  data  base  is  created  from  element 
data  and  special  elements  (springs,  dashpots,  etc)  are  inserted  at  appropriate 


points  in  the  data  base.  All  element  information  (shell  sections,  beam 
sections,  distributed  loads,  etc.)  are  interpreted  and  resolved  onto  the 
element  data  base  and  the  short  list  data  base  (volume  4).  Also,  the  nodes 
are  assigned  with  internal  numbers  based  on  the  element  ordering. 

Step  5.  Build  node  point  data  base  -  Node  information,  such  as  point  loads, 
boundary  conditions,  transformations,  etc.,  are  interpreted  and  resolved  onto 
the  node  point  data  base  and  the  short  list  data  base  (volume  5). 

Step  6.  Build  step-  dependent  part  of  short  list  data  base  -  the  step- 
dependent  part  of  the  short  list  data  base  (volumes  6  -  9)  is  created  for  each 
analysis  step. 

Step  7.  Space  allocations  for  the  MAIN  -Based  on  the  global  sizing  values  for 
the  MAIN  analysis  program,  space  allocations  for  the  MAIN  are  made.  The 
parameters  for  the  array  pointers  (starting  addresses  of  arrays)  are  stored  in 
common  blocks.  The  maximum  sizing  is  limited  by  a  parameter  MCOREM,  which  is 
currently  set  to  be  45,000  words.  For  problems  exceeding  this  capacity,  a 
new  value  may  be  assigned  to  MCOREM  in  the  main  driver  of  PRE. 

Step  8.  Plot  undeformed  mesh. 

Step  9.  Build  communication  data  base  -  In  order  to  pass  all  necessary 
information  (values  of  array  pointers,  control  variables,  etc.)  from  the  PRE 
onto  the  MAIN,  a  communication  data  base  is  created  and  copied  on  disk  or  tape 
unit  #23.  It  consists  of  common  blocks,  element  data  base,  node  point  data 
base,  material  point  data  base,  and  short  list  data  base. 

Finally,  a  driver  for  the  MAIN  is  written  onto  a  disk  or  tape  unit  #28.  Two 
typical  drivers  are  shown  as  the  following: 
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a)  For  Stress  Analysis 

DOUBLE  PRECISION  SINT 
COMMON  /  CSP  /  SINT  (45000) 

CALL  EXEC* 

STOP 

END 

♦Executive  module  for  stress  analysis 

b)  For  Heat  Transfer  Analysis 

DOUBLE  PRECISION 

COMMON  /  CSP  /  SINT  (45000) 

CALL  HEXEC  A* 

STOP 

END 

*  Executive  module  for  heat  transfer  analysis. 

The  driver  is  compiled  into  an  object  module  and  kept  on  a  temporary 
file.  After  linking  the  driver  to  the  object  module  of  the  MAIN,  the  analysis 
of  MAIN  Program  starts. 

3.2  MAIN  Program. 

The  MAIN  contains  a  number  of  modules  which  perforin  the  actual 
engineering  analysis.  In  the  current  version  of  ABAQUS,  there  are  three  major 
analysis  branches: 

•  Stress/strain  analysis 

•  Offshore  analysis  (not  active  in  reviewer's  version) 

•  Heat  transfer  analysis 

•  Eigenvalue  analysis 
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For  a  typical  analysis  branch,  several  modules  are  executed  in  a  loop 
many  times  over  due  to  the  incremental  (or  time-history)  solution  approach 
being  used.  The  main  loop  consists  of  three  prime  calculation  stages;  i) 
material  calculations,  ii)  element  calculations,  and  iii)  solution  of 
matrix  equations.  Detailed  descriptions  of  the  modules  in  the  MAIN  are  given 
in  ABAQUS  system  manual  [9].  In  this  section,  we  shall  review  ABAQUS' 
organization  by  considering  the  following; 

•  Input  Step 

•  Procedure  Library 

•  Material  Library 

•  Element  Library 

•  Equation  Solver 

•  Arithmetic  Library 

3.2.1  Input  Step 

In  accordance  with  the  type  of  analysis  requested,  an  executive 
manager  is  called  to  read  the  data  on  communication  data  base  from  disk  #23 
sequentially.  The  buffer  size  of  communication  data  base  was  set  to  be  512  in 
word  length.  All  data  in  the  buffer  area  of  communication  data  base  are  then 
copied  to  appropriate  common  blocks  and  data  bases.  Tape  unit  for  each  data 
base  is  also  defined  as  shown  in  Table  3.1.  The  core  storage  arrangement  is 
problem  dependent.  For  a  relatively  small  problem,  the  buffer  area  is  divided 
into  several  segments  and  a  work  space.  Each  segment  is  occupied  exclusively 
by  a  data  base.  Un  the  other  hand,  for  large  problems  several  data  bases  may 
share  a  coimnon  segment  in  the  buffer.  All  the  pointers  for  the  buffer 
division  are  saved  in  comon  blocks.  For  example,  in  the  material  generation 


3-6 


F 


Table  3.1  Tape  Arrangement  in  the  MAIN  Program 


Tape 

Number 

Sequential/ 

Random 

Description 

1 

S 

Eigenvalue  data  base  /  half-step  residual 

file 

2 

S 

Decomposed  equation  file 

3 

S 

Element  nickname  data  base  /  Quasi  Newton 

data 

base* 

4 

S 

Temporary  file  for  node  point  data  base 

5 

S 

Read  input 

6 

S 

Write  input 

8 

S 

FILE  output  data  base 

12 

S 

Restart  data  base  * 

19 

s 

Element  operator  data  base  * 

20 

s 

Element  operator  data  base 

21 

s 

Material  point  data  base  * 

22 

s 

Material  point  data  base 

23 

s 

Communication  data  base 

24 

s 

Material  stiffness  data  base  /  temporary  element 
Output  file  /  solution  nickname  data  base/post 
plot  file 

25 
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Element  data  base 
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29 
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Node  point  data  base 
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Element  mass  matrix  data  base 

*0ne  for  read 

1  -  one  for 

write 
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module  pointers  are  stored  in  the  common  block  called  "CMATGI",  whereas  in  the 
element  generation  module  pointers  are  in  ''CELRI“.  Most  of  the  array  pointers 
^  are  defined  in  the  PRE  program  although  some  of  the  pointers  are  calculated  in 

the  MAIN  whenever  they  are  needed.  In  view  of  this,  two  immediate  comments 
can  be  made:  i)  If  any  array  allocations  are  to  be  redefined  in  the  MAIN, 

^  related  changes  on  the  definitions  of  the  pointers  must  also  be  made  in  the 

PRE  program,  and  ii)  changes  become  difficult  if  the  data  structure  of  the 
data  bases  is  involved. 

3.2.2  Procedure  Library 

Currently,  there  are  four  executive  managers  which  control  the  procedure 

I 

options,  namely 

•  Static  Analysis  (EXECS) 

•  Dynamic  Analysis  (EXECD) 

•  Natural  Frequency  Calculations  (EXECEV) 

•  Heat  Transfer  Analysis  (HEXEC) 

The  relationship  of  the  procedure  library  for  stress  analysis,  with  the  input 
module  in  the  MAIN  is  shown  in  Fig.  3.1.  Macro  flow  for  each  executive 
manager  is  almost  the  same.  Typical  steps  for  static  stress  analysis  are 
outlined  as  the  following  (Fig.  3.2): 

i)  Initialize  step  and  set-up  increment  information, 

ii)  Calculate  material  stiffness  D  (i.e.  stress-strain  relations) 
for  each  material  point, 

iii)  Calculate  element  stiffness 

iv)  Solve  the  system  of  equations  by  the  frontal  method 
;  v)  Update  the  material  stiffness  for  the  case  of  nonlinear 


material s 


Fig.  3.1.  Modules  of  Procedure  Library 


Set  up  Increment 
Information 
(Inc) 


Calculate  Material 
Stiffness 
(MATGEN) 


Calculate  Element 
Stiffness 
,  EIGEN) 


Solve  System 
Equations  by 
Frontal  Method 
(Solve) 


Update  Node  Point  Data 
Base  (2u),  Calculate 
Strain,  Write  it  to 
material  DB. 


Update  Material  Data 
Base 
(MATGEN) 


Calculate  Stresses 
and  Residuals 
(ELRES) 


/Solution 

^Convergence 


Update  Nodal 
and  Element 
Variables 


Print  Output 
and  Nodal 
Output 


r^Step 
Fini shed 


(.Next  increment) 


Fig.  3.2  Macro  Flow  of  Static  Stress  Analysis 
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vi)  Calculate  stresses  and  residual  forces  at  integration  points 

vii)  If  the  solution  has  converged,  update  the  nodal  point  data 
base 

viii)  Print  output  for  element  and  nodal  information 

ix)  Proceed  the  analysis  for  the  next  time  increment. 

3.2.3  Material  Library 

The  library  is  controlled  by  the  module  driver  called  f'lATGEN.  The 
pointers  for  defining  the  workspace  are  stored  in  the  common  block  'CMATGI'. 
The  material  stiffness  data  base  pointers  and  other  related  information  to  the 
data  bases  are  initialized  in  the  subroutine  MINIT.  For  every  material  point 
(which  may  differ  from  the  element  integration  point),  material  type  is 
checked  and  the  corresponding  control  variables  are  defined.  Six  different 
types  of  material  models  (or  routines)  can  be  called  to  generate  the  material 
stiffness  matrix  D; 

•  MATELA  -  Linearly  elastic  materials  with  isotropic,  orthotropic  and 

anisotropic  properties. 

•  MATEXP  -  Linearly  elastic  thermal  expansions  with  isotropic, 

orthotropic  and  anisotropic  properties. 

•  MATELG  -  Linearly  elastic  materials  for  general  sections  (options 

available  for  beam  and  shell  elements). 

•  MATEXG  -  Linearly  elastic  thermal  expansions  for  general  sections. 

•  MATCRP  -  Creep  models  (strain-hardening,  time-hardening  or  user  input 

creep  law)  and  irradiation  induced  swelling. 

•  MATPLA  -  Plasticity  model  (von  Mises  plasticity  with  isotropic  or 

kinematic  hardening  law).  The  effects  or  rate-  and 
temperature-dependence  core  also  available. 


3-11 


Also,  other  material  options,  such  as  density,  conductivity,  specific  heat, 
etc.  for  heat  transfer  analysis,  are  defined  in  module  HEXEC. 

Two  stages  of  calculations  are  involved  for  material  generation; 
a)  Material  stiffness  calculations,  and  b)  Material  data  updating.  For  the 
material  stiffness  calculations,  the  following  steps  are  executed  for  all 
material  points: 

i)  Set  up  material  data 

ii)  Calculate  the  thermal  strain,  if  it  is  required, 

iii)  Calculate  initial  strain  for  creep  and  swelling,  whenever 
appropriate 

iv)  Compute  tangent  modulus. 

v)  Write  the  data  on  the  material  stiffness  data  base. 

For  the  updating  stage,  executions  include: 

i)  Compute  stress  increments. 

ii)  Insert  new  state  on  material  point  data  base  record, 

iii)  Write  the  data  on  the  material  point  data  base. 

A  flow  chart  for  the  material  calculations  is  shown  in  Fig.  3.3.  Both  the 
material  stiffness  and  material  point  data  bases  provide  necessary 
communication  information  between  the  material  library  (module  MATCEN)  and  the 
element  library  (module  ELGEN).  It  is  seen  that  in  ABAQUS  the  material 
library  is  completely  independent,  in  terms  of  coding  organization,  from  the 
element  library.  That  is,  the  material  models  are  uniformly  available  to  all 
element  types.  This  feature  is  different  from  ADINA  [11]  in  which  the 
material  library  is  a  sub-module  of  each  element  type  [2]. 


! 


08HK:  iVaterial  Stiffness  Data  Base. 

SELUP:  Update  Special  Element  Information. 


Fig.  3.3  Flow  Chart  for  Material  Generations 
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3.2.4  Element  Library 

This  library  is  commanded  by  the  routine  ELGEN  (module  driver).  Its 
primary  purpose  is  to  calculate  element  stiffness  and  other  element  data  and 
then  update  the  record  on  the  element  stiffness  data  base.  Calculations  are 
made  in  a  loop  for  all  elements.  The  major  calculation  steps  involve  the 
fol lowing : 

i)  Call  ELSET  to  read  element  information 

ii)  Kead  data  for  all  material  points  of  the  element  from  material 

stiffness  data  base  to  workspace. 

iii)  Initialize  element  stiffness  matrix  [K]  and  element  residual  vector 

lRe>- 

iv)  Form  element  strain  -  nodal  displacement  transformation  matrix  [B], 

v)  For  dynamic  analysis: 

a.  Calculate  the  mass  matrix  [M]  for  the  first  time  increment 

b.  Add  the  effect  of  initial  conditions  to  the  right  hand  side  force 
vector  {Rg}. 

c.  Form  effective  stiffness  matrix  [S]. 

vi)  Check  the  value  of  control  variable  NREF  for  appropriate  adjustment 

of  stiffness  matrix  [K]  or  [S],  and  force  vector  {Rg}  : 

a.  Additional  terms  due  to  distributed  loads  or  foundation  stiffness. 

b.  Coordinate  transformation  from  local  to  gobal  system. 

c.  The  effect  of  boundary  condition. 

d.  Equations  for  multi-point  constraints. 
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vii)  For  the  first  iteration  of  an  increment,  write  {Rg}  to  node  point 
data  base.  Otherwise,  read  {Rg>  from  node  point  data  base, 
viii)  Write  element  stiffness  data  base  record  consisting  of  [S]  and  [R]. 

3.2.5  Arithmetic  Library 

This  library  consists  of  several  utility  packages  which  are  accessible  to 
all  the  modules  in  ABAQUS.  These  packages  are: 

•  Matrix  addition 

•  Matrix  subtraction 

•  Matrix  multiplication 

•  Element  stiffness  calculation;  /  [B]^[D]  [B]. 

V 

•  Calculation  of  element  residual  force  vector, 

{Rgl  =  {Re>  -  /  {a}  dv 

V 

•  Copy  a  double-precisioned  real  array  to  an  integer  array. 

•  Copy  an  integer  array  to  a  double-precisioned  real  array. 

•  Cholesky  decomposition  of  a  matrix. 

•  Q-R  algorithm 

•  Vector  normalization,  etc. 

Obviously,  by  consolidating  all  the  utility  subroutines  in  a  library,  ABAQUS 
has  avoided  any  repetitive  coding  as  it  was  commonly  done  in  other  programs. 

Finally,  to  show  the  inter-relationship  of  various  modules  in  the  MAIN 
program,  a  chart  is  drawn  in  Fig.  3.4.  The  hierarchical  level  of  the  modules 
is  that  any  module  can  only  call  other  modules  inside  of  its  circle.  One  may 


note  that  the  interaction  between  the  element  library  and  the  material  library 
is  made  through  material  stiffness  and  material  point  data  bases. 

Commentary : 

1)  Separation  of  the  entire  program  into  PRE  and  MAIN  provides  apparent 
advantages ; 

a)  Short  turn-around  time  -  Running  the  PRE  program  only  takes  a  few 
seconds  or  a  few  minutes  of  CPU  time,  which  has  high  priority  in 
execution. 

b)  Smaller  core  requirement  -  In  PRE,  all  real  variables  are  in  single 
precision  and  it  requires  no  more  tl.an  750  K  bytes  core  size  to 
analyze  a  fairly  sizable  finite  element  problem. 

The  PRE  can  be  run  either  separately  or  consecutively  with  the  MAIN. 

2)  The  code  was  very  modilarly  constructed  for  ease  of  extension  or 
modification  by  the  developers. 

3)  The  codings  of  ABAQUS  are  generally  difficult  to  follow  by  a  non-developer 
since  there  are  very  few  comment  cards  in  the  program  for  explanation. 

4)  The  material  library  is  an  independent  module  from  the  element  library. 
Therefore  any  existing  material  models  in  ABACUS  are  uniformly  available 
to  all  element  types. 

5)  Although  ABAQUS  was  modularly  constructed,  to  add  a  new  material  model  or 
element  type  into  the  program  by  a  user  is  by  no  means  obvious.  Not  only 
one  has  to  add  subroutines  corresponding  to  the  new  material  model  or  new 
element,  but  also  changes  have  to  be  made  in  the  root  subroutines  of  the 
main.  Moreover,  a  number  of  subroutines  in  the  PRE  must  be  modified. 
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3.3  Data  Base  Design 

ABAQUS  utilizes  the  data  base  concept  in  both  PRE  and  MAIN  programs.  In 
this  concept,  the  data  are  separated  into  several  files,  each  file  is  handled 
independently.  In  the  program  a  master  array  called  INTS  is  allocated  and  it 
is  divided  into  data  base  page  pools  (or  buffer  areas)  and  working-space.  All 
data  files  are  stored  in  the  page  pools  with  pre-ass igned  sizes.  The  sizing 
decision  is  made  in  the  PRE  program.  For  example,  the  data  base  structure 
used  in  the  PRE  is  shown  in  Fig.  3.5.  The  lengths  of  the  page  pools  will  not 
remain  fixed,  they  can  be  redefined  during  the  analysis.  All  page  pools  are 
allocated  in  the  fast  core  area.  If,  for  a  particular  data  file,  the  page 
pool  is  not  large  enough,  secondary  storage  device  (either  sequential  or 
random  access  file)  is  utilized. 

Access  of  each  data  base  is  controlled  by  its  manager.  That  is,  whenever 
information  needs  to  be  calculated  or  updated,  the  data  base  manager  will 
transfer  the  data  from  page  pool  area  to  working-space.  After  the 
calculations  are  completed,  data  are  then  copied  back  to  the  page  pool  area. 
Each  page  pool  is  a  physical  record  and  it  is  further  separated  into  a  number 
of  fields.  The  field  is  the  smallest  unit  that  can  be  accessed  by  the  data 
base  manager  in  ABAQUS.  Moreover,  each  field  consists  of  several  pieces  of 
information.  The  process  of  transferring  data  between  the  page  pool  area  and 
work-space  is  shown  in  Fig.  3.6. 

Clearly,  the  data  base  manager  performs  the  following  important  functions: 

i)  Transfer  data  from  the  data  base  (page  pool)  to  work  space 

ii)  After  the  data  has  been  processed,  copy  it  back  to  the  page  pool  area, 

iii)  If  secondary  storage  device  is  used,  perform  necessary  I/O  operation, 
then  transfer  data  as  i  or  ii. 


3-18 


Array  INTS 

Restart  data  base  (read/write) 

Restart  data  base  (write/read) 

Node  list  file,  history  data  base 

Element  list  file.  Material  point  data  base,  etc. 

Node  set  file.  Communication  data  base. 

Element  set  file 

Input  data  base 
Element  data  base 
Short  list  data  base 

Output  file  data  base 

Node  point  data  base 

Expanded  node  point  data  base 
Input  data  base  manager 

Working  space 
Notes : 

^  NGTRC  is  defined  as  1024. 

*  Temporary  working  area  for  the  arrays  indicated. 


Fig.  3.5.  Buffer  Areas  of  Data  Bases  in  Pre-Program 


3-20 


From  the  above  discussion,  we  may  summarize  the  following  characteristics: 

i)  By  using  the  data  base  design  concept  and  separate  work-space  for 
processing  data,  all  data  are  well  protected  from  accidental 
over-writing  or  deletion. 

ii)  By  keeping  independent  data  bases,  program  modularity  can  be  easily 
achieved.  For  example,  the  material  library  is  an  independent  module 
from  the  element  module;  interactions  between  the  two  modules  are  made 
through  the  material  stiffness  and  material  point  data  bases, 

iii)  Data  integrity  is  maintained  throughout  the  analysis;  that  is,  data 
bases  contain  only  processed  data,  all  calculations  or  updatings  are 
performed  in  the  v/orkspace. 

iv)  The  usage  of  secondary  storage  devices  Is  kept  minimum,  it  varies  with 
the  size  of  a  pro)lem. 

Nevertheless,  there  a”e  two  drawbacks  which  may  warrant  future 
improvements : 

i)  Pre-defined  fixed  sizes  of  data  bases  may  cause  unnecessary  wasting  of 
storage,  which  will,  in  turn,  increase  I/O  swaps.  By  doing  this,  some 

of  the  data  bases  (in  core)  may  have  been  exhausted,  whereas  others 

I 

still  have  excess  space  or  are  not  being  used  at  all. 

ii)  By  using  separate  workspace  for  data  processing,  large  amounts  of  data 
copying  are  involved  and  this  defintely  increases  the  CPU  time 
significantly. 


4.  FUNCTIONAL  UESCKIPTION 

In  this  section,  an  overview  of  the  analysis  capabi lities*of  ABAQUS  and  a 
description  of  solution  methodology  are  given.  It  is  noted  again  that  all 
discussions  are  centered  on  the  Version  3  of  ABAQUS. 

4.1  Analysis  Procedures 

ABAQUS  is  designed  principally  for  the  nonlinear  static  and  dynamic 
analysis  of  structures,  and  nonlinear  steady  and  transient  analysis  of  heat 
transfer  or  conduction  problems.  Linear  analysis  for  the  corresponding 
problems  can  also  be  done  as  a  special  case  of  the  nonlinear  analysis.  The 
Version  3  of  ABAQUS  has  the  following  three  basic  procedures: 

•  Static  stress  analysis 

•  Dynamic  analysis 

•  Heat  transfer  analysis  (both  steady  state  and  transient). 

•  Eigenvalue  Analysis 

4.1.1  Static  Stress  Analysis 

ABAQUS  has  fairly  extensive  nonlinear  analysis  capabilities.  For  stress 
analysis,  three  major  sources  of  nonlinearity  are  included; 

i)  Material  nonlinearity  -  Possible  material  models  are  to  be  discussed 
later  in  the  material  library. 

ii)  Geometric  nonlinearity  -  Only  the  effect  of  large  rotation  but  small 
strain  is  considered  in  the  current  version, 

iii)  Boundary  nonlinearity  -  The  user  may  impose  gap  or  no-gap,  and 

friction  boundary  conditions  between  two  bodies.  Classical  Coulomb 
friction  is  assumed. 

The  user  may  request  any  combinations  of  the  aforementioned 
noni inearities. 

♦Additional  analysis  capabilities  of  Version  4  -  ABAQUS  are  listed  in  Appendix  C. 
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The  solution  method  employed  in  ABAQUS  is  an  incremental  approach  with 
tangent  stiffness.  For  each  loading  increment,  (full)  Newton's  method  is  used 
to  solve  nonlinear  equilibrium  equations.  The  obvious  advantage  of  Newton's 
method,  over  the  modified  Newton  or  quasi -Newton ,  is  that  it  has  better 
convergence  rate  and  consequently  the  method  can  handle  a  wider  range  of 
nonlinear  problems.  On  the  other  hand,  this  procedure  is  rather  expensive 
since  the  structural  stiffness  matrix  (or  heat  conduction  matrix  for  heat 
transfer  analysis)  has  to  be  reformed  and  triangulari zed  for  every  iteration 
cycle. 

In  ABAQUS,  the  history  of  a  problem  is  divided  into  three  stages;  step, 
increment,  and  iteration  cycle.  A  problem  may  consist  of  one  or  several 
steps;  typically  each  step  is  simply  a  change  from  one  load  to  another  (say 
from  static  to  dynamic  load).  A  step  then  consists  of  a  number  of  increments. 
Within  each  increment,  several  iteration  cycles  are  involved.  Loading  history 
and  solution  convergence  are  controlled  by  the  following  input  parameters: 

INC  -  Maximum  number  of  increments  in  a  step 

CYCLE  -  Maximum  number  of  iteration  cycles  allowed  in  an  increment 

PTOL  -  Force  tolerance  for  checking  solution  convergence  during 
iterations. 

MTUL  -  Moment  tolerance  for  checking  solution  convergence  during 
iterations. 

The  user  manual  recommends  the  values  of  PTOL  and  MTOL  being  from  1  to  10%  of 
typical  applied  forces  and  reaction  forces  in  a  structure.  The  reason  for 
using  such  physical  quantities  for  convergence  control  is  based  on  the  as¬ 
sumption  that  the  engineer  should  know  well  enough  about  his  model  and  be  able 
to  specify  appropriate  values  for  PTOL  and  MTOL.  Contrary  to  this  opinion, 
the  reviewer  feels  that  it  is  difficult  to  define  a  realistic  PTOL  or  MTOL 
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for  d  complicated  structure  since  the  analyst  often  does  not  have  any  notion 
as  to  the  order  of  reaction  forces  or  moments  a  priori. 

For  loading  increment  control,  the  user  has  two  options  to  choose  from: 
i)  Direct  step  control  -  user  specifies  the  loading  increments  . 
ii)  Automatic  step  -  The  program  will  choose  its  own  incrementation 
scheme  based  on  specified  tolerance,  i.e.  PTUL  and  MTOL. 

As  explained  in  the  manual,  direct  step  control  is  suggested  in  cases  where 
the  problem  behavior  is  well  understood.  Otherwise,  automatic  step  is 
preferred,  especially  for  start-up  of  a  new  nonlinear  problem.  Generally, 
automatic  step  requires  much  more  CPU  time  than  direct  step  for  a  given 
problem. 

To  exercise  the  automatic  step,  two  additional  input  parameters  are 
requi red: 

NUMBER  -  Suggested  number  of  increments  of  uniform  size  needed  to 
complete  the  solution. 

CUTMAX  -  The  maximum  number  of  times  current  increment  size  may  be 
subdivided. 

The  program  will  proceed  the  analysis  based  on  a  suggested  time  increment  At 
(total  time)/NUMBER,  and  monitor  the  maximum  force  and  moment  residuals. 
Starting  at  iteration  cycle  #4,  solution  convergence  is  projected  at  n-th 
iteration  that 

n  =  i  +  itn  (PT0L/R(i))/4n{R(i  (4-1) 

where 

i  =  i-th  itertion  cycle 

R(^j  =  Maximum  force  or  moment  residual  corresponding 
to  i-th  iteration 
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Once  n  is  calculated,  two  possibilities  exist: 

a)  n  £  CYCLE,  convergence  is  expected  and  hence  continue  iterations. 

b)  n  >  CYCLE,  the  suggested  time  increment  is  too  large  and  cut-back  of 
increment  size  is  made,  i.e. 

New  At  =  0.25  (Old  at) 

The  solution  is  then  started  over  again  from  the  end  of  the  previous 
increment. 

On  the  other  hand,  if,  in  two  successive  increments,  convergence  is 
achieved  in  less  than  CYCLE/2  iterations,  the  increment  size  is  increased  by 
25%  for  the  next  increment,  with  a  minimum  number  of  increments,  NUMBER. 

According  to  the  benchmark  problems  run  by  the  reviewer  (to  be  shown  in 
Sections  5  and  6),  the  automatic  step  control  worked  quite  well  in  terms  of 
getting  a  solution  successfully.  However,  in  most  cases  this  procedure  is 
extremely  time-consuming  and  it  is  impractical  to  solve  a  large  scale 
nonlinear  problem  using  automatic  step  exclusively  from  the  beginning  to  the 
end. 

4.1.2  Dynamic  Analysis 

With  the  inclusion  of  inertia  effect,  the  equilibrium  equations 
corresponding  to  finite  element  approximation  can  be  expressed  as 

Mir  +  !-  ?■  =  0  (4-2) 

where 

M  =  Consistent  mass  matrix 

u  =  d^u/dt^;  u,  nodal  displacement  vector 
T  =  Internal  force  vector 


P  =  External  force  vector 
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For  nonlinear  problems,  internal  force  vector  I  is  related  to  the  tangent 
stiffness  matrix  of  the  finite  element  idealization  which  may  include  the 
effects  of  geometric  as  well  as  material  nonlinearities. 

To  solve  Eq.  (4-2)  numerically,  two  stages  of  approximations  are 
involved;  Newton's  method  of  iterations  as  employed  in  static  analysis  and 
direct  integration  in  time.  ABAQUS  adopted  the  a  -  method  due  to  Hilber  and 
Hughes  [12,  13]  which  is  an  implicit  integration  scheme.  In  this  method,  the 
equilibrium  equations  are  written  as 

^  ^t+At  *  *  “)(^t+At  ■  ^t+At)  -  a  (Tt  -  ^t)  “  0  (^-3) 

In  addition,  Newmark's  formulae  are  used  to  approximate  the  displacement  and 
velocity  vectors  at  time  t+At, 

\+At  =  \  ^  ^t  ^  ^  ®VAt^ 

•  • 

IJt+At  =  Ut  +  At  ((1  -  y)  Ut  +  T  (4-5) 

witn 

1 

6  =  4  (1-  a)2 
1 

Y  -  Y  "  ® 

-  ^  C  a  < 

The  main  characteristics  of  the  a  -  operator  are; 


X 
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i)  It  is  unconditionally  stable  in  linear  dynamic  analysis 

ii)  It  offers  controllable  algorithmic  damping  by  specifying  the  a-values 

iii)  The  operator  is  more  accurate  in  the  lower  modes,  as  compared  to  the 
Wilson  6  -  method  [14].  For  large  At/T  ratio  (where  T  =  period),  it 
ensures  strong  dissipation  in  the  higher  modes, 

iv)  For  the  dynamic  effect,  a  =  -  0.05  was  recommended. 

Similar  to  the  static  analysis,  the  user  has  the  options  to  proceed  the 
time  history  analysis  by  using  either  direct  step  control  or  automatic  step 
control.  For  direct  step,  solution  convergence  is  determined  by  comparing  the 
residual  force  and  moment  at  t+At  with  PTUL  and  MTOL.  On  the  other  hand,  for 
automatic  stepping  convergence  is  determined  by  checking  the  residual  force  at 
the  half-step  between  t  and  t+At,  or  called  half-step  residual  Hi/2* 
quantity  is  calculated  as  follows. 

The  accelerations  u  are  assumed  to  vary  linearly  over  the  time 
interval  [t,  t+At] 

Ut  “  (  I  -  T  )  Uj.  +  T  ,  0  <  X  <  1  (4~6) 

Based  on  the  above  approximation,  by  direct  integration  one  can  find  the 

•  ** 

expressions  for  Ux.  and  u^  in  terms  of  the  corresponding  quantities  at  t 
and  t+At  .  Then  from  Eq.  (4-2),  the  residual  force  vector  at  any  t  within  the 
step  is  defined  by 

Tr(T)  =  TT  u  (t  )  +  r  (t  )  -  F  (t)  (4-7) 

For  convergence  check,  a  half-step  residual  chosen  as  a  measure 

against  an  input  tolerance  called  HAFTOL.  The  value  Rifz  ’s  the  magnitude  of 
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the  laryest  entry  (absolute  value)  in  the  vector  R  for  x  =  1/2.  The  basis 
for  such  a  choice  is  purely  empirical. 

To  specify  an  appropriate  tolerance  value  for  HAFTOL  can  be  a  challenging 
task  to  the  user.  According  to  the  manual,  the  following  values  for 
HAFTUL  are  recommended:  If  P  is  a  typical  magnitude  of  real  forces  in  a 
problem,  then 

i)  For  high  accuracy,  HAFTUL  *  lO"^  P 

ii)  For  moderate  accuracy,  HAFTOL  *  10“^  P 

iii)  For  coarse  increment  size,  HAFTOL  '  P 
As  will  be  seen  in  Section  6,  for  a  given  problem  the  magnitude  of  ^\/2 
is  rather  sensitive  to  the  a-  value.  If  a  =  0,  very  large  half-step  residual 
was  obtained  from  the  analysis  and  consequently,  this  results  in  severe 
cut-back  of  increment  size  in  the  automatic  time  stepping.  Therefore,  when 
the  automatic  step  control  is  optioned,  a  =  0  should  not  be  used  for  the  most 
cases.  Even  with  a  =  -0.05  or  other  nonzero  value,  some  significant 
oscillations  of  ^\/2  observed  at  the  initial  stage  of  the  dynamic 
response.  In  summary,  the  following  commentary  can  be  made: 

i)  Automatic  time  stepping  is  an  extremely  useful  feature  for  nonlinear 
dynamic  analysis,  particularly  for  start-up  of  an  analysis.  However, 
the  automatic  stepping  should  not  be  used  exclusively  for  the  entire 
his.  ry  of  a  problem,  since  it  usually  leads  to  more  CPU  time.  Some 
degree  of  user's  intervention  or  judgment  is  probably  necessary  if 
one  is  to  keep  the  analysis  time  within  a  reasonable  limit, 
ii)  The  magnitude  of  half-step  residual  quite  sensitive  to  the  a- 

value.  One  should  avoid  the  usage  of  o  =  0  if  automatic  time  step  is 
exercised.  As  recommended  in  ABAUUS  user  manual,  a  =  -0.05  is 
probably  the  best  damping  ratio  for  practical  purposes. 
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4.1.3  Heat  Transfer  Analysis 

The  third  general  procedure  in  Version  3  of  ABAQUS  is  the  heat  transfer 
analysis.  It  is  primarily  for  solid  body  heat  conduction  with  temperature 
dependent  conductivity,  internal  energy  (including  latent  heat  effects),  and 
convection  and  radiation  boundary  conditions.  Thermo-mechanical  coupling  was 
ignored.  The  heat  transfer  process  can  be  either  steady  state  or  transient. 

For  time  integration  of  transient  heat  transfer  analysis,  a  modified 
Crank-Nicholsen  operator  (backward  difference  algorithm)  was  adopted  in 
conjunction  with  Newton's  method  of  iterations. 

In  terms  of  coding  organization,  the  heat  transfer  part  is  an  independent 
module,  parallel  to  the  stress  analysis  module.  This  module  has  its  own 
element  library  and  material  library.  The  output  of  heat  transfer  analysis 
are  the  nodal  temperatures,  which  can  be  directly  used  for  thermal  stress 
analysis. 

4.2  Element  Library 

In  ABAQUS,  there  are  two  major  categories  of  elements:  Elements  for 
stress  analysis  and  elements  for  heat  transfer  analysis. 

Stress  Analysis: 

•  Truss  elements 

•  2/3  Solid  elements 

•  3/D  Solid  elements 

•  Beams 

•  Shell  elements 


Elbows 
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Heat  Tra.isfer  Analysis: 

•  Link  elements 

•  Plane  and  axisymmetric  solids 

•  3/D  Solid  elements 

A  brief  summary  o'‘  the  above  elements  is  given  in  Table  4.1.  Each  of  the 
element  types  is  outlined  as  below. 

4.2.1  Truss  Elements 

These  are  one-dimensional  bar  elements  which  can  have  axial  deformation 
only.  Two  different  truss  elements  are  available:  2-node  truss  with  linear 
displacement  and  3-node  truss  with  quadratic  displacement.  In  3/D  space,  each 
node  has  three  translational  degrees  of  freedom,  u^,  Uy,  u^. 

4.2.2  2/D  Solid  Elements 

These  are  first  order  (4-node  quadrilateral)  and  second  order  (8-node 
quadrilateral)  isoparametric  elements.  Moreover,  the  elements  are  divided 
’  'to  three  sub-classes  in  accordance  with  the  state  of  deformations:  a) 
plane  stress,  b)  plane  strc»in,  and  c)  axisymmetric  deformation.  For  each 
element,  the  orders  of  numerical  integration  for  evaluating  element  stiffness, 
stress/strai n  output,  consistent  mass  matrix,  etc.,  are  fixed  in  the  program. 
Reduced  integration  scheme  is  available  for  8-node  elements  which  are  useful 
for  such  cases  where  near  incompressible  materials  or  significant  plastic 
deformations  are  encountered. 

4.2.3  3/D  Solid  Elements 

These  are  equivalent  to  2/D  solid  elements  except  that  they  are  defined 
in  three-dimensional  space.  The  elements  include  8-node  brick  with  linear 
displacement,  20-node  brick  with  quadratic  displacement,  and  20-node  brick 
with  quadratic  displacement  and  reduced  integration. 
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Table  4.1  Element  Library  Summary 


No.  of  Nodes 

OOF  per 

No.  of  Integration  Points 
Stiffness  and  Mass  and 

Element  Tvoe 

per  Element 

node 

Stress/Strain  Forces  * 

Truss: 


Linear  Displacement  2 

3 

1 

2 

Quadratic  Displace.  3 

3 

2 

3 

2/D  Solid: 

(Plane  stress,  plane 
strain,  axi symmetric) 

Bilinear  Displacement  4 

2 

2x2 

2x2 

Biquadratic  Displace.  8 

2 

3x3 

3x3 

Biquadratic  with 

Reduced  Integration  8 

2 

2x2 

3x3 

3/0  Solid: 

Linear  Displacement  8 

Quadratic  Displac.  20 

Quadratic  with 
Reduced  Integration  20 


3  2x2x2  2x2x3 

3  3x3x3  3x3x3 


3 


2x2x2  3x3x3 


Beam: 


Planar  Beam 

Straight  Linear  2 
Curved  Quadratic  3 
Curved  Cubic  2 


3 

3 

3 


1+ 

2 

3 


2 

3 

3 


Space  Beam: 

Straight  Linear  2 
Quadratic  Linear  3 
Cubic-Curved  2 
Cubic-Straight  2 


6 

6 

6 

6 


2 

3 

3 


2 

3 

3 

3 


Shell  Elementst: 

Linear  Displacement 
Quadratic  Displace. 


4 


6 


I 


2x2 
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Table  4.1  Element  Library  Summary  (continued) 


Element  Ty£e 


No.  of  Integration  Points 
No.  of  Nodes  DOE  per  Stiffness  ana  Mass  and 
per  Element  node  Stress/S^rain  Forces  * 


Elbow: 


Linear  Displacement  2 

Linear  with 

Ovalization  2 

Quadratic  Displace.  3 

Quadratic  with 

Ovalization  3 


6 


1^ 


6 

6 


1 

2 


6  2 


2 

2 

3 

3 


Heat  Transfer  Elements 
Link  Elements: 


2- Node  Link  2 

3- Node  Link  3 


Specific  Heat 
Conductivity  and  Flux _ 


1  2  2 

1  3  3 


Plane  and  Axisymetric 

4-Node  Element  4 

8-Node  Element  8 


1  2x2  2x2 

1  3x3  3x3 


3/0  Elements 

8-Node  Element  8 

20-Node  Element  20 


1  2x2x2  2x2x2 

1  3x3x3  3x3x3 


Notes : 

*  Integrations  including  consistent  mass  matrix,  body  forces,  specific  heat 
and  flux,  surface  pressure  or  flux,  etc. 

+  For  planar  beam,  space  beam  and  elbow  elements,  the  indicated  integration 
points  are  along  the  axis  of  the  element.  Integr^^’on  scheme  over  the  cross 
section  of  the  beam  or  elbow  is  specified  separate./  in  the  user  manual. 

f  For  shell  elements,  reduced  integration  scheme  is  mandatorily  enforced  in 
the  program. 
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4.2.4  Beams 

Several  beam  elements  with  a  wide  variety  of  cross-sectional  shapes  are 
available  in  ABAQUS.  Cross  sections  include 

•  General  section  -  User  defines  section  properties 

•  Pipe  section 

•  Box  section 

•  Circular  solid  section 

•  I  -  section 

•  Rectanyular  solid  section 

•  Hexagonal  box  section 

The  beam  elements  are  further  divided  into  two  different  types  according  to 
their  spacial  configurations. 

A.  Planar  Beams 

•  2  node  straight  beam  -  linear  interpolation  function  and  elastic 
transverse  shear  deformation 

•  3  node  curved  beam  -  quadratic  interpolation  function  and  elastic 
transverse  shear  deformation 

•  2  node  curved  beam  -  cubic  interpolation  function  and  elementary 
beam  theory 

B.  Space  Beams 

In  addition  to  those  equivalent  to  the  planar  beams,  the  fourth  entry  is 


2  node  straight  beam  -  cubic  interpolation  function  and  elementary 
beam  theory 


4-13 


Numerical  integration  for  evaluating  element  properties  is  performed  somewhat 
differently  from  other  element  types:  Along  the  beam  axis,  Gauss  integration 
scheme  is  used.  Over  the  cross  section,  however,  default  integration  is  the 
Newton-Cotes  method  whereas  Gauss  integration  is  optional  for  some  of  the 
cross  sections.  The  Newton-Cotes  method  is  preferred  for  maximum  (or  minimum) 
bending  stresses  occurred  near  the  outer  fibers  of  the  beam. 

4.2.5  Shell  Elements 

ABAOUS  has  fairly  sophisticated  shell  elements  for  conducting  three- 
dimensional  thin  shell  analysis.  Two  different  elements  of  the  same  type  are 
available  for  use:  4-node  and  8-node  elements;  each  node  has  6  -  degrees  of 
freedom,  i.e.  3  -  translations  and  3  -  rotations  referring  to  a  global 
Cartesian  coordinates.  The  shell  equations  were  derived  from  a 
two-dimensional  classical  shell  theory  [15],  differing  from  the  so-called 
degenerated  3/D  continuum  theory  [16]. 

The  shell  is  defined  by  its  middle  surface  which,  in  turn,  is 
approximated  by  the  nodal  coordinates  and  isoparametric  functions.  All 
displacement  components  and  rotations,  and  the  geometry  of  the  middle  surface 
are  interpolated  by  the  same  'serendipity'  functions.  The  direct  strain  in 
the  thickness  direction  of  the  shell  is  assumed  to  be  zero.  Since  the  element 
is  Specifically  intended  for  thin  shell  analysis,  Kirchhoff  assumption  with  no 
transverse  shear  deformation  was  imposed;  i.e.  normal  vectors  to  the  middle 
surface  of  the  shell  always  remain  normal  during  deformation.  This  assumption 
is  enforced  by  using  a  penalty  function  at  a  set  of  (uniformly)  reduced 
integration  points. 

Based  on  its  formulations  and  sample  problems  tested,  it  appears  that  the 
shell  elements  in  ABAQUS  have  the  following  characteristics: 


i)  The  elements  are  very  reliable  for  solving  a  wide  range  of  thin  shell 
structures. 

ii)  It  gave  better  numerical  behavior  for  shells  of  small  thickness,  as 
compared  to  the  serendipity  element  with  reduced  integration, 

iii)  The  shell  must  have  uniform  thickness  within  the  element. 

4.?. 6  Elbow  Elements 

Elbows  are  often  used  in  the  design  of  piping  systems  for  nuclear  power 
plants,  petro-chemical  facilities,  etc.  It  is  known  that  elbows  achieve  their 
flexibility  through  a  shell-like  behavior,  resisting  bending  actions  with 
significant  cross-sectional  oval izations.  For  elastic  elbows,  siffness  matrix 
based  on  the  3/D  straight  beam  theory  is  used  in  conjunction  with  flexibility 
factors  and  the  stresses  are  corrected  by  using  appropriate  intensity  factors 
L17].  However,  for  inelastic  elbows  no  such  simple  approach  is  available 
and  one  has  to  result  in  other  numerical  means.  In  this  connection,  a  number 
of  elbow  elements  have  been  proposed  [18-21]  and  they  generally  suffer  two 
drawbacks: 

i)  Element  formulations  were  made  on  the  over-simplified  theory  and 
consequently  its  applications  become  restrictive, 
ii)  Computational  cost  is  very  high. 

From  its  theoretical  considerations,  the  elbow  elements  in  ABAgUS 
probably  represent  one  of  the  most  attractive  types  for  nonlinear  (material) 
analysis.  The  element  formulations  were  derived  on  the  basis  of  the  following 
assumptions: 
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•  Only  small  deformations  are  considered. 

•  Uvalization  of  the  pipe  cross  section  is  allowed  and  the  effect  of 
warping  may  also  be  included. 

•  The  total  deformation  of  the  pipe  section  is  considered  to  be  the 
linear  superposition  of  strains  caused  by  beam  actions  and  strains  by 
ovalization  end  uniform  radial  expansion  of  the  pipe  section. 

•  Interpolation  functions  for  beam  deformation  are  similar  to  the 
elementary  beam  elements. 

•  Interpolation  for  cross-section  deformation  in  the  axial  direction 
can  be  either  linear  or  quadratic  polynomial  function.  Up  to  six 
Fourier  modes  (sine  and  cosine  functions)  can  be  specified  for 
deformation  associated  with  ovalization  of  the  pipe  section. 

4.2.6  Heat  Transfer  Elements 

As  mentioned  previously,  only  solid  body  is  considered  for  heat  transfer 
analysis  in  ABAQUS.  This  is  sufficiently  general  for  the  geometric 
description  of  most  application  problems.  Therefore,  the  elements  include: 

•  l/D  link  (bar)  elements  -  2  nodes  or  3  nodes 

•  2/D  plane  and  axi symmetric  elements  -  4  nodes  or  8  nodes 

•  3/0  solid  elements  -  8  nodes  or  20  nodes 

4.3  Material  Library 

In  the  present  version  of  ABAQUS,  analysis  is  limited  to  large  rotations 
but  small  strains  for  which  case  the  second  Piol a-Ki rchhof f  stress  and  Green's 
strain  are  used  for  stress  and  strain  measures,  respectively.  In  this  way, 
any  rotational  effects  of  constitutive  properties  due  to  the  deformation  of 
material  coordinates  are  already  accounted  for. 
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Basically,  there  are  four  different  types  of  material  ttwdels  for  stress 
analysis  in  conjunction  with  heat  transfer  models: 

Stress  Analysis 

a)  Elastic  models  -  Isotropic,  orthotropic  or  fully  anisotropic 
materials 

b)  von  Mises  Plasticity  -  Isotropic  or  Hill's  anisotropic  materials 
with  plastic  incompressibility 

-  Isotropic  strain  -  hardening 

-  Kinematic  strain  -  hardening 

c)  Creep  models  -  Isotropic  or  Hills'  anisotropic  materials, 
uniaxial  creep  equation  : 


A  (/* 


where 


c  =  creep  strain  rate 
a  =  equivalent  stress 


t  =  time 


and  A,  n,  m  are  material  constants  which  may  be  functions  of 
temperature.  Or,  user  may  specify  his  own  creep  equation  through 
user  subroutine.  Creep  hardening  options  include 
Time  hardening  law 
Strain  hardening  law 

d)  Volumetric  swelling 

This  model  was  prepared  primarily  for  the  stress  analysis  of  fuel 
elements  in  reactor  core  where  irradiation  induced  strain  (swelling) 
occurs  in  flux  fields.  The  rate  of  swelling  is  defined  as  a  function  of 
temperature  and  other  field  variables,  such  as  time  and/or  radiation 


flux.  The  user  has  the  option  to  input  the  swelling  data  in  tabular  form 
or  specify  the  swelling  law  by  a  user  subroutine. 


It  is  noted  that  all  material  constants  involved  in  the  above  models  can 
vary  with  temperature.  Moreover,  thermal  stresses  can  be  analyzed  through  the 
option  'EXPANSION'. 

Heat  Transfer  Analysis 

a)  Conductivity  -  Isotropic,  orthotropic  or  fully  anisotropic 

b)  Specific  heat  -  This  option  is  needed  for  transient  heat  transfer 
analysis.  Specific  heat  is  defined  as  a  function  of  temperature. 

As  it  has  been  discussed  before,  each  material  model  is  uniformly  available  to 
all  element  types  in  ABAQUS. 

4.4  Mesh  Generation 

ABAQUS  has  a  simple  version  mesh  generation,  i.e.,  generating  a  mesh 
of  uniform  grid.  The  mesh  generation  option  is  provided  in  nodal  point  data 
and  element  data  input,  led  by  two  keyword  cards: 

*  NGEN  -  Incremental  node  generation 

*  ELGEN  -  Incremental  element  generation 

The  program  can  generate  11  the  intermediate  nodes  along  a  straight  line  or 
curved  line  (parabolic  or  circular)  by  specifying  the  nodal  numbers  and 
coordinates  of  two  end  nodes.  If  a  parabola  is  chosen,  the  user  must  supply 
the  coordinates  of  the  mid-point  on  the  curve  between  the  two  end  points.  On 
the  other  hand,  if  a  circular  arc  is  used,  the  user  must  define  the 
coordinates  of  the  center  of  the  circle  or  the  number  of  center  as  one  node. 

In  addition,  the  option  ELGEN  is  used  to  generate  a  series  of  elements  with 
uniform  element  and  nodal  number  increments.  For  this  purpose,  first  the  user 
has  to  specify  a  master  element  and  associated  nodal  numbers. 
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Then  the  program  will  use  the  master  element  to  generate  successive  elements 
of  the  same  type  with  constant  element  and  nodal  increments. 

Therefore,  it  is  clear  that  the  version  3-ABAQUS  does  not  yet  have  any 
sophisicated  mesh  generation  features  beyond  the  uniform  grid. 

4.5  Kinematic  Conditions 

In  addition  to  the  fixed  or  prescribed  displacement  boundary  conditions, 
four  kinematic  constraining  features  are  available  in  AQAQUS: 

•  EQUATION  -  User  defined  multi -point  constraint  by  a  linear  equation. 

•  FRICTION  -  Friction  between  node  pairs. 

•  GAP  -  Possible  contact  and  gap  conditions  between  node  pairs. 

•  MPC  -  General  multi-point  constraints. 

Each  option  is  defined  by  a  leading  keyword  card  with  additional  input  data. 
These  features  are  briefly  discussed  as  follows. 

4.5.1  EQUATION 

This  feature  provides  an  option  to  specify  multi-point  constraints  in  a 
linear  equation  form 

A  u  +  A  u„  +  ...  +  A  U  =  0  (4-8) 

*  ‘  ^  ^  n  n 

where  A^,  A^,  ...,  Ap  are  user  defined  values;  u^,  u^,  ....  Up  are  nodal 
degrees  of  freedom,  which  are  specified  by  node  numbers  and  their  degrees  of 
freedom.  This  option  is  an  altenative  to  the  MPC  option  with  more  flexibility 
in  the  manner  of  constraining  the  nodal  degrees  of  freedom. 

4.5.2  FRICTION 

This  is  an  option  to  define  the  interaction  effect,  namely  friction, 
between  the  contacting  surfaces  of  two  bodies  and  it  is  usually  used  in 


4-19 


conjunction  with  GAP  option.  The  friction  action  is  represented  by  the 
Coulomb  model  for  which  the  friction  property  is  a  function  of  the  normal 
force  acting  at  the  interface  of  contact.  For  this  option,  the  user  must 
supply  the  following  data 

•  The  pair  of  node  number  in  contact 

•  Coefficient  of  friction 

•  Stiffness  in  Stick  Condition 

•  Closure  distance  (if  there  is  a  gap) 

•  Normal  force  (if  gap  is  always  closed). 

•  Direction  cosines  of  closure  distance 

•  Direction  cosines  of  first  friction  direction. 

Although  the  program  considers  the  friction  effect,  it  does  not  print  out  slip 
deformation  between  two  nodes.  Moreover,  the  friction  relationship  is  limited 
to  small  deformation. 

4.5.3  GAP 

This  condition  allows  pairs  of  nodes  to  be  in  contact  or  separated;  it  is 
equivalent  to  a  bar  -  contact  -  element.  The  user  specifies  the  direction  and 
closure  d  of  two  nodes,  which  represents  the  initial  gap  between  two  bodies. 
Then  the  program  monitors  the  relative  displacement  of  the  two  gap  nodes. 

When  the  relative  displacement  reaches  the  value  d,  the  gap  is  closed  and 
subsequently  a  constraint  is  imposed. 

The  program  handles  both  static  and  dynamic  contact  problems.  For 
dynamic  contact,  a  set  of  impact  equations  are  solved  to  obtain  the  velocity 
and  acceleration  jumps  of  those  nodes  in  contact  during  the  instant  before  and 
after  the  contact  has  occurred.  Then,  the  initial  velocities  and  acceleration 
of  the  contacting  nodes  can  be  calculated  and  the  usual  time  stepping 
integration  proceeds. 


I 

i 
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4,5.4  MPC  -  Multi  Point  Constraints 

This  option  is  useful  for  the  transition  region  between  a  fine  mesh  and 
coarse  mesh  of  a  finite  element  model  layout.  As  seen  in  Fig.  4.1,  the  mesh 
is  refined  from  coarse  to  fine  with  two  different  cases:  i)  linear 
displacement  elements  ii)  quadratic  displacement  elements.  In  either  case, 
the  node  P  on  ab  line  or  nodes  P^  and  on  abc  lines  must  be  constrained  to 
have  campatible  displacement  field  for  the  elements  on  both  side  of  the  node. 
If  these  nodes  are  not  constrained,  significant  solution  errors  could  occur, 
especially  when  plasticity  is  spread  in  the  transition  zone  [22].  In  ABAQUS, 
the  MPC  option  is  available  for  both  2/0  and  3/0  solid  elements 

4.6  User  Subroutines 

In  addition  to  the  standard  data  input,  certain  data  can  be  specified 
through  user  subroutines.  This  option  provides  added  flexibility  and 
convenience  for  the  program  to  generate  data  which  follow  specialized  equation 
or  patterns.  The  user  subroutines  available  are: 

CREEP  -  To  calculate  creep  rate  from  a  specific  creep  law. 

DFLUX  -  To  define  a  non-uniform  distributed  flux  in  heat  transfer 
analysis. 

DISP  -  To  specify  the  magnitude  of  displacement  for  all  degrees  of 

freedom  of  the  nodes  where  the  boundary  conditions  are  defined. 

DLOAD  -  To  define  non-uniform  distributed  loads  in  an  element. 

FILM  -  To  define  non-uniform  film  coefficient  and  associated  sink 
temperatures  for  heat  transfer  analysis. 

MPC  -  To  define  a  multi-point  constraint  similar  to  but  more  general 
than  the  option  in  Section  4.5.4. 
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SWELL  -  To  define  a  time  dependent  volumetric  swelling  law. 

Although  the  usage  of  some  of  the  subroutines  is  fairly  obvious,  specific 
examples  for  each  option  should  be  given  in  the  manual  to  povide  better 
clarity  on  the  parameters  involved  in  the  subroutines. 

4.7  Input  and  Output 

Data  input  in  ABAQUS  is  organized  in  two  major  sections: 

Model  input  -  specifying  geometric  definitions  of  finite  element  model 

such  as:  nodal  number,  coordinates  of  nodes,  element-node 
connectivity,  constraint  conditions,  etc.;  material  model 
information;  and  plotting  of  undeformed  geometry. 

Hisotry  Input-These  include  data  related  to  analysis  procedure,  loading 
history,  output  control  and  solutin  control  (e.g.  direct  or 
automatic  stepping  method),  etc. 

Input  data  are  divided  into  distinct  sections  and  each  section  has  a  leading 
keyword  card.  A  keyword  card  begins  with  a  notation  (e.g.  *  ELEMENT), 
also  contains  optional  parameters.  Following  the  keyword  card  are  the  data 
cards.  On  any  data  card,  free  or  fixed  format  may  be  used.  It  is  the 
reviewer’s  opinion  that  ABAQUS  input  has  the  following  merits: 

i)  Input  sequent  was  logically  layed  out  and  therefore  it  is  easy  for  the 
user  to  follow  the  input  manual. 

ii)  Separation  of  data  sections  by  keyword  cards  enhances  the  user's  ability 
to  scan  input  data  quickly  and  spot  possible  errors. 

iii)  Free  format  input  is  most  convenient  to  use. 
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With  respect  to  the  output  of  ABAQUS,  four  types  of  output  are  available: 
•  Printed  output 
.  Plotting 

•  File  output  for  post-processing  analysis 
.  Restart  file  output 

Our  discussion  is  focused  on  the  printed  output. 

In  ABAQUS,  the  printed  output  is  divided  into  two  sections;  namely,  in 
the  PRE  and  the  MAIN. 

In  PRE: 

•  Input  echo  (or  input  card  image) 

•  A  summary  of  various  options  processed  by  the  program 

.  Processed  input  data  -  element  definitions,  material  description, 
node  definitions,  analysis  step  information,  loads,  mesh  plotting, 
etc. 

•  Work  space  allocation 
In  MAIN: 

•  Element  stresses  and  strains 

«  Total  nodal  displacements,  velocities  and  accelerations  (for  the 
case  of  dynamic  analysis). 

•  Nodal  residual  or  reaction  forces 

.  Information  about  iteration  and  solution  convergence 
The  user  has  the  option  to  suppress  some  of  the  printout  through  the  input 
section  '*PRINT'.  Among  various  options,  for  example,  the  user  may  choose: 
i)  Printout  can  be  requested  at  every  k-th  interval. 
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ii)  The  user  has  the  option  to  print  out  the  element  data  for  a  specific 
number  of  elements  which  are  designated  in  one  ELSET. 

iii)  The  printout  of  element  variables,  such  as  stress  components, 
stress  variants,  total  strain  components,  creep  strains,  etc, 
can  be  specified. 

iv)  Printout  of  node  set,  total  displacements,  velocities,  accelera¬ 
tions,  etc.  may  also  be  optioned. 

It  is  apparent  that  ABAQUS  provides  a  very  flexible  printout  option 
to  the  user,  which  represents  one  of  the  friendly  features  of  the  code. 
Nevertheless,  two  negative  points  about  its  printed  output  are  mentioned 
for  consideration  of  future  improvement; 

i)  Printout  of  real  numbers  are  in  P-format  which  is  the  mixture 
of  F-  and  E-  format.  Although  the  use  of  such  format  is  well 
intentioned,  it  does  not  have  a  good  appearance  for  report  pur¬ 
pose.  The  reviewer  still  prefers  the  use  of  E-  format  for  printing 
real  numbers. 

ii)  Stress  and  strain  components  are  not  defined  in  the  printout. 


Consequently,  it  is  not  convenient  for  the  user  to  identify  their 
meanings  unless  the  input  manual  is  referenced. 
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5.  Verification  Problems 

Eleven  verification  problems  were  run  to  check  the  execution  status  of 
ABAQUS.  The  purpose  of  such  verification  exercises  is  purely  to  ensure  the 
program,  which  was  developed  on  a  CDC  Cybernet  System,  is  functioning  properly 
on  the  IBM  370-158  computer  with  double  precision  arithmetic  on  floating  point 
numbers.  The  verification  problems  represent  different  types  of 
nonlinearities,  but  are  of  routine  nature  in  that  they  do  not  present  any 
numerical  difficulty  in  obtaining  convergent  solutions.  Each  problem  is 
presented  in  a  standardized  format  including: 


i) 

Problem  statement 

ii) 

Ueomet  ry 

iii) 

Nonlinearity 

iv) 

Material  properties 

V) 

Loading 

Vi) 

Finite  element  irodel 

vii ) 

Analysis  procedure 

viii ) 

Results 

For  later  discussion,  several  parameters  being  used  in  ABAQUS  are  defined  as 
fhe  following: 

R  =  Maximum  reaction  force 
E  =  Young's  modulus 
V  =  Poisson's  ratio 

Oq  =  Initial  yield  stress 

Ep  =  Plastic  modulus 

p  =  Density 
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5.1  Elastic-plastic  Deformation  of  a  Cantilever  Beam 

Geometry:  A  short  cantilever  beam  with  a  rectangular  cross  section. 

Length  =  91.4  cm.  Height  =  30.5  cm.  Width  =  2.54  cm. 

Nonlinearity:  Material  nonlinearity  with  small  deformation 
Material  properties:  Elastic-perfectly  plastic  material 
E  =  68.95  mpa,  v  =0 

Oq  =  0*59  mpa 

Loading:  A  concentrated  transverse  load  P  applied  at  the  free  end. 

Model:  Two  different  models  were  used. 

A.  lO-beam  (B22)  elements 

B.  20-eight  node  plane  stress  (CPS8)  elements. 

Procedure:  Automatic  load  incrementing 

Results:  The  load  deflection  response  at  the  free  end  is  plotted  in  Fig. 

5.1  for  both  models  A  and  B  in  conjunction  with  an  independent  solution  by 
Felippa  [23],  It  is  seen  that  the  elastic  responses  of  model  A  and  model  B 
agree  with  that  of  Felippa.  However,  the  results  obtained  from  beam  elements 
show  some  discrepancy  near  the  transition  zone  between  elastic  and  plastic 
deformations.  This  is  probably  due  to  insufficient  integration  order  being 
used  in  the  code  for  this  problem.  Nevertheless,  the  limit  loads  from  iTOdels 
A  and  B,  and  reference  [23]  are  about  the  same. 


5.2  A  Cantilever  Beam  Subjected  to  An  End  Moment. 

Geometry:  A  slender  cantilever  beam  with  a  rectangular  cross  section 


(Fig.  5.2) 

Length  =  10m,  Height  =  0.1m,  Width  =  1  m. 

Nonlinearity:  Large  rotation  and  linear  material 

Material  Properties:  Linearly  elastic  material,  E  =  1.2  x  10^  KN/m^ 

Loading:  A  moment  applied  at  the  free  end 

Model:  8-beam  (B21)  elements 

Procedure:  Automatic  load  incrementing 
Results: 

After  few  trial  runs,  it  was  found  that  a  model  consisting  of  8  beam 
elements  is  sufficient  to  represent  the  bending  action  of  the  beam,  including 
large  deformation  response.  The  moment  was  applied  in  such  a  manner  that  the 
beam  was  deformed  from  its  straight  shape  to  a  complete  circle.  Three 
different  sets  of  solution  parameters  were  specified  to  determine  the  solution 
sensitivity: 


Case 

PTOL 

MTOL 

NUM 

1 

0.3  lb 

0.3(in-lb) 

10 

2 

0.3 

0.3 

30 

3 

0.03 

0.3 

10 

4 

0.1 

0.1 

10 

The  moment  vs.  vertical  and  horizontal  deflections  are  shown  in  Figs.  5.2 
and  5.3,  respectively.  All  results  correlate  with  those  obtained  by  Horrigmoe 
[24]  quite  well  except  the  Case  1  in  which  some  localized  numerical  instability 
was  indicated  for  the  horizontal  deflection  (see  dotted  line  in  Fig.  5.3). 
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Fig.  5.3).  This  instability  was  removed  by  either  reducing  the  force 
tolerance  or  increasing  the  irtaximum  number  of  iterations.  Fig.  5.4  shows  the 
deflection  curves  of  the  beam  corresponding  to  different  values  of  applied 
moment . 

5.3  Static  Analysis  of  A  Pipe  Contact  Problem 
Geometry:  (see  Fig.  5.5) 

Pipe:  =  300  in.  L2  =  100  in. 

L3  =  50  in.  *-4  =  120  in. 

R=  45  in.  r  =  15in.  t  =  l^in. 

Restraint  L5  =  27  in  Area  =  26  in^  . 

’Gap  h  =  3  in. 

Nonlinearity:  Nonlinear  boundary  condition  (gap  and  contact)  and  elastic- 
plastic  material 

Material  Properties: 

Pipe  -  E  =  27  X  103  ^si  oo  =  28  ksi 

Restraints  -  E  =  30  x  103  ksi  oq  =  38  ksi 

Loading:  A  vertical  concentrated  force  P  =  -50K  was  applied  at  the  free  end. 

Model:  The  model  consists  of  (see  Fig.  5.6) 

Pipe  -  7  straight  beam  (B21)  elements 
1  curved  beam  (B22)  element 
All  elements  have  pipe  cross  sections. 

Restraints  -  2  truss  (C102)  elements. 

Procedure:  Automatic  load  incrementing. 

This  problem  was  selected  from  the  ABAQUS  example  problem  manual  [10]  to 
test  the  analysis  capability  of  gap  elements.  Three  different  cases  of 
solution  parameters  were  specified: 
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Case 

PTOL 

NUM 

1 

O.IK  (.2%)* 

20 

2 

O.IK  (.2X) 

2 

3 

50K  (100%) 

5 

*  Percentage  of  the  total  applied  load  or  maximum  reaction  force. 

The  change  of  gap  size  and  the  development  of  contact  force  in  the  supports 
were  plotted  against  the  applied  load  shown  in  Figs.  5.7  and  5.8, 
respectively.  It  is  seen  that  the  results  for  all  cases  are  almost  identical. 
The  only  difference  was  the  CPU  time,  i.e.  Case  1=2  min. -32  sec.  Cases  2  &  3 
10  sec.  For  this  problem,  the  source  of  nonlinearity  is  primarily  due  to  the 
contact  action  between  the  pipe  and  supports.  Ideally,  two  load  increments 
would  be  sufficient  to  complete  the  analysis;  namely,  load  up  to  the  point  of 
initial  contact  and  the  remaining  total  load.  In  all  cases,  the  residual 
force  calculated  for  each  loading  incement  was  very  small,  well  below  the 
smallest  tolerance  specified. 

5.4  An  Elastic-Plastic  Spherical  Cap 

Geometry:  A  spherical  shell  with  initial  imperfection  (axi symmetric) 
as  shown  in  Fig.  5.9. 

Two  different  cases  of  dimensions  were  considered. 

Case  A:  R  =  0.8177  in.  Rj  =  1  .1432  in.  di  =  30<’ 

h  =  0.0104  in.  a  =  0.297  in. 

Case  B:  R  =  0.825  in.  Rj  =  1.1511  in.  di  =  20® 

h  =  0.025  in.  a  =  0.267  in. 

The  cap  was  fixed  along  its  edge. 


Gap  (in.) 
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Load  P  (%) 

Fig.  5.7.  Change  of  Gap  Size  vs.  Applied  Load 


Contact  Force  (Kips) 


5-i3 


—  Case  1 
•  Case  2 
X  Case  3 


Load  P  (%) 

Fig.  5.8.  Development  of  Contact  Force  at  Support  #1  vs.  Applied  Load 
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Nonlinearity:  Elastic-plastic  material  and  sii'ill  djfornation. 

Material  properties:  Bi-linear  elastic-plastic  material 
E  =  1.8  X  103  ksi  V  =  0.3 

oo  =  78  ksi  Ep  =  1.1  X  103  ski 

Loading:  Uniform  normal  pressure  P 

Model:  10  -  axisymmetric  2/0  eight-node  (CAX8R)  elements  with  reduced 

integration. 

Procedure:  Automatic  load  inci'ementing 
PTOL  =  10  lb.,  NUM  =  5 

Results: 

This  problem  was  previously  analyzed  by  Kao  [25]  using  the  finite 
difference  method  for  axisymmetric  shells.  The  cap  has  an  initial 
imperfection  which  is  also  axisymmetric.  The  load  vs.  deflections  at  the 
crown  for  Case  A  and  Case  B  are  shwon  in  Figs.  5.10  and  5.11,  respectively. 

In  general,  ABAQUS  results  agree  with  those  in  [25],  except  that  a  small 
discrepancy  is  seen  for  Case  B. 

5.5  Elastic-Plastic  Deformation  of  A  Plate 

Geometry:  A  square  plate  with  a  span  of  40  in  and  uniform  thickness  of  0.4  in. 
Two  different  support  conditions  were  considered: 

A.  Clamped  along  all  edges 

B.  Simply  supported 

Nonlinearity:  Elastic-plastic  material  and  small  deformation 
Material  properties:  Perfectly  elastic-plastic  material 
E  =  30  X  103  Ksi  V  =  0.3 

oq  =  30  ksi 


0.0104  In 


0.825  In 
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Loading;  Uniform  pressure  p 

Model:  A  4  x  4  mesh  consisting  of  16  shell  (S8R)  elements  with  reduced 

i ntegration 

Procedure:  Both  automatic  and  direct  load  incrementing  procedures  were  used. 
Results : 

For  the  clamped  boundary  conditions,  the  non-dimensional ized  load  vs. 
maximum  deflection  of  the  plate  obtained  from  ABAQUS  was  plotted  in  Fig.  4.12 
as  compared  with  existing  solutions.  ABAQUS  result  agrees  closely  with  that 
of  Horrigmoe  and  Eidsheim  [26],  and  the  ultimate  load  of  the  plate  coincides 
with  upper  bound  obtained  by  Hodge  and  Belytschko  [27]. 

For  the  simply  supported  plate,  ABAQUS  solution  correlates  with  that  in 
[26]  quite  well  until  the  load  reaches  the  ultimate  value.  Some  numerical 
disparity  is  clearly  seen  in  Fig.  5.13.  It  is  noted  that  the  analysis  was 
performed  using  automatic  load  incrementing  with  the  solution  parameters: 

PTQL  =  25  lb,  (2%  of  the  maximum  reaction  force),  NUM  =  10.  Initially,  ABAQUS 
solution  followed  that  of  Horrigmoe  and  Eidsheim  quite  closely  all  the  way 
near  the  ultimate  load.  Then,  convergence  became  difficult  and  ABAQUS  cut 
back  the  load  increments  in  order  to  achieve  solution  convergence.  As  a 
result,  gradual  stiffening  of  the  plate  can  be  seen  in  Fig.  5.13  until  the 
load  approaches  the  upper  bound  solution,  analysis  was  interrupted  due  to 
numerical  instability  of  the  plate  (or  reaching  ultimate  load).  Similar 
stiffening  behavior  occurred  even  when  PTUL  was  reduced  with  greater  iteration 
allowance.  It  is  believed  that  the  stiffening  effect  was  caused  by  some 
unknown  numerical  difficulty  in  ABAQUS. 


PLVM. 


t  »  thickness  0.4  in 
*  central  deflection 
P  *  uniform  load 

D  »  Et»/12Ci.^a) 


^  Yield  line 
theory24,0 


24.9 

.Upper  and  lower 
bound.  Hodge  and 
Belytschko  [27] 


initial  yielding  1S.Q 


Horrigmoe  and  Eidsheim  [26] 


ABACUS 


0.05 


0. 10 


0.15 


“c  ■ 


0.20 


Fig.  5.13,  Load-Oeflectfon  Curve  of  an  Elastic-Plastic 
Simply  Supported  Plate 
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5.6  A  Circular  Plate  with  Plasticity 

Geometry:  A  clamped  circular  plate  with  dimensions 

Radius  a  =  10  in.  ,  thickness  t  =  1  in. 


Nonlinearity:  Elastic-plastic  material  and  small  deformation 

Material  Properties:  A  Ramberg-Osgood  stress-strain  curve  was  assumed  to 

represent  the  strain-hardening  behavior  of  the  material. 


where 


o 


e  =  Total  effective  strain 

a  =  Total  effective  stress 

o  =  24  ksi 
n  =  6.67 


Also  E  =»  lO**  ksi,  V  =  0.33  ,  Oq  =  16  ksi 
Loading:  Uniform  pressure  p 

Model:  From  symmetry,  one  quarter  of  the  plate  was  modeled  by  12-  shell 

(S3R)  elements  as  shown  in  Fig.  5.14. 

Procedure:  Direct  load  incrementing. 

Results: 

Several  runs  were  made  by  specifying  different  solution  tolerances,  i.e. 
PTOL  and  MTOL,  ranging  from  0.7  to  20%  of  the  maximum  reaction  forces.  It  was 
found  that  the  solution  corresponding  to  the  maximum  tolerance  differs  from 
that  of  the  minimum  tolerance  by  about  8%  at  the  highest  loading  level.  On 
the  other  hand,  when  PTOL  and  MTOL  £  5%  of  the  maximum  reaction  force,  the 
solutions  can  hardly  be  distinguishable. 
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Results  obtained  from  ABAQUS  are  compared  with  those  from  Kao  [25]  and 
Popov  et  al  [28]  in  Fig.  5.5  and  good  agreement  can  be  seen. 

5.7  Large  Deformation  of  a  Cylindrical  Shell 

Geometry:  A  cylindrical  shell  (shown  in  Fig.  5.16)  has  the  dimensions 
R  =  2540  mm  i  =  2540  mm 

h  =  3.175  fnn  0  =  0.1  rad. 

All  edges  are  clamped. 

Nonlinearity:  Large  rotation  and  linearly  elastic  material 

Material  Properties:  E  =  3.1  KN/mm^,  v  =  0.3 

Loading:  Uniform  pressure  normal  to  the  surface  of  the  shell. 

Model:  A  4  x  4  mesh  consisting  of  16  shell  (S8R)  elements. 

Procedure:  Direct  loads  incrementing,  PTOL  =  ION  (3%  R),  MTOL  =  25  N.M., 

Results: 

The  purpose  of  this  problem  is  to  test  the  large  deformation  feature  of 
the  shell  element  in  ABAQUS,  including  both  the  pre-  and  post-buckling 
responses.  This  shell  has  been  a  rather  popular  choice  by  several 
investigators  [24,  29  and  30]  due  to  its  stable  behavior  during  the  buckling 
stage  so  that  no  numerical  difficulty  was  encountered.  As  seen  in  Fig.  5.16, 
the  post-buckling  response  shows  some  noticeable  differences  among 
various  investigators.  One  possible  explanation  is  the  inclusion  of 
deformation  -  dependent  pressure,  which  was  considered  in  our  ABAQUS  run. 


O  V  Popov  et  al .  [28] 

Fig.  5.15.  Deflection  Curves  of  a  Circular  Plate 
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024  6  8  10  12 

Central  deflection,  W^(nin) 


a 

A 


} 


Horrigmoe,  Different  Elements  [24] 


-  Dhatt  [30] 

Sabir  &  Lock  [29] 
+  ABAQUS 


Fig.  5.16.  Load-Deflection  of  a  Cylindrical  Shell. 
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5.8  Post-Buckling  of  a  Cylindrical  Shell  Subjected  to  a  Point  Load 

Geometry:  The  dimensions  are 

R  =  2540  mm  1  =  2540  mri 

h  =  12.7  rnm  6  =  0.1  rad 

The  two  longitudinal  edges  were  hinged  and  immovable,  whereas 
the  curved  ends  were  free. 

Nonlinearity:  Large  rotation  and  linearly  elastic  material. 

Material  Properties:  E  =  3.1  KN/mm^  ,  v  =  0.3 

Loading:  A  concentrated  vertical  force  applied  at  the  center. 

Model:  A  4  x  4  mesh  consisting  of  16-shell  {S8R)  elements. 

Procedure:  Automatic  load  incrementing 
Results: 

Although  the  geometry  of  this  problem  is  very  similar  to  the  preceding 
benchmark,  its  post-buckling  behavior  is  quite  different  due  to  the  simply 
supported  boundaries  and  concentrated  load.  In  order  to  carry  on  the  analysis 
beyond  the  buckling  point  of  the  shell,  a  displacement  boundary  condition, 
instead  of  load  control,  must  be  specified  at  the  center  of  the  shell.  The 
load  vs.  deflection  response  is  shown  in  Fig.  5.17.  The  solution  from  ABAQUS 
is  almost  identical  to  that  of  Horrigmoe  [24]. 

5.9  Dynamic  Response  of  a  Cantilever  Beam 

Geometry:  A  cantilever  beam  with  a  rectangular  cross  section 
length  =  10  in.,  height  =  1  in.,  width  =  1  in. 

Nonlinearity:  Two  cases  were  analyzed. 

A.  Small  deformation 

B.  Large  deformation 


Central  load,  P(KN) 
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o 


o 
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Material  properties:  Linearly  elastic  material  was  assumed  for  both  cases, 
i.e.  E  =  1.2  X  lO**  psi,  u  =  0.2 
p  =  10“^  Ib-sec^/in**. 

Loading:  The  beam  is  subjected  to  a  step-loading  uniform  pressure  of  2.85 
psi,  the  pressure  was  equally  distributed  on  the  top  and  bottom 
surfaces  of  the  beam  as  shown  in  Fig.  5.18. 

Model:  The  cantilever  was  modeled  by  5  beam  (B21)  elements. 

Procedure:  Two  different  time  increments  were  analyzed,  namely  Atj  = 

4.5  X  10~^  sec  and  At2  =  9  x  10“5  sec. 

Results: 

The  purpose  of  this  problem  was  to  test  the  dynamic  algorithms  of  ABAQUS. 
The  same  problem  was  previously  analyzed  by  ADINA  [11].  The  maximum 
deflection  histories  for  small  and  large  deformation  cases  are  shown  in  Figs. 
5.19  and  5.20,  respectively.  It  is  seen  that  ABAQUS  results  agree  with  those 
of  ADINA. 


aunssajj 


Normal ized  Time 


Fig.  5.20.  Deflection  History  of  a  Cantilever  Beam  With 

Large  Deformation 
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6.  ADVANCED  EVALUATION 

Several  advanced  benchmark  problems  were  selected  to  test  the  numerical 
characteristics  of  ABAQUS.  These  include 

i)  Selection  of  solution  parameters,  such  as  PTOL  or  HAFTOL,  etc. 

ii)  Sensitivity  of  solution  to  convergence  control, 

iii)  Analysis  time  vs.  procedure  option  and  the  range  of  convergence 
tolerance  specified. 

iv)  Execution  time  of  ABAQUS  vs.  ADINA. 

The  benchmarks  selected  are  those  with  simple  geometries  so  that  data 
preparation  effort  was  kept  at  a  minimum.  These  problems  are  relatively 
complex  in  terms  of  their  numerical  natures. 

6.1  An  Elastica 

Elastica  is  a  well-known,  classical  problem,  for  which  the  closed-form 
solution  is  available  in  reference  [31].  Several  investigators  have  taken 
this  problem  as  a  benchmark  to  verify  the  large  deformation  algorithm  of  their 
finite  element  programs  [24,  32].  Although  successful  finite  element 
solutions  were  obtained  by  previous  authors,  caution  must  be  given  to  this 
problem  from  the  numerical  standpoint. 

A  bar  of  rectangular  cross-section  is  subjected  to  uniaxial  load  as 
shown  in  Fig.  6.1.  Material  was  assumed  to  be  linearly  elastic;  Young's 

I 

modulus  E  =  30,000  ksi ,  Poisson's  ratio  v  »  0.3.  A  beam  element  B22  was 
chosen  to  represent  the  bending  action  of  the  bar  and  the  corresponding  finite 
element  model  is  given  in  Fig.  6.1.  In  order  to  activate  the  buckling  of  the 
bar,  an  initial  imperfection  in  the  y-ordinates  of  the  nodes  was  introduced 
according  to  y  =  6  •  cos  (n  x/2L),  6  »  0.1. 
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Several  computer  runs  were  made  using  ABAQUS  and  It  was  found  that  the 
solution  of  this  problem  was  rather  sensitive  to  three  parameters  tested: 
a)  mesh  size,  b)  increment  size,  and  c)  convergence  tolerance  for 
iterations. 

Initially,  a  relatively  coarse  mesh  consisting  of  5  B22  beam  elements 
was  used  to  tnodel  the  beam.  This  model  appeared  to  be  satisfactory  when 
compared  with  elementary  beam  theory  within  the  small  deformation  range.  In 
running  ABAQUS,  a  prescribed  nodal  displacement  of  150"  at  the  right  end  of 
the  bar  and  automatic  incrementing  procedure  were  specified.  The  deflection 
shapes  for  several  runs  by  varying  the  number  of  solution  increments  (NUM)  and 
convergence  tolerance  (i.e.  PTOL  and  MTOL)  are  tabulated  in  Table  6.1.  It  is 
seen  that  when  PT0L=  MTOL  =  1,  and  NUM  =  1,  pure  compressive  mode  of 
deformation  was  obtained  from  the  analysis.  While  for  PTOL  =  MTOL  =  10,  and 
NUM  =  5  and  20,  respectively,  the  corresponding  deformation  mode  was  switched 
to  a  direction  opposite  to  that  of  the  induced  initial  imperfection.  Whereas, 
for  other  cases  the  results  were  close  to  the  analytical  solution  for  the 
nodal  displacement  greater  than  100",  but  considerable  differences  were  found 
for  smaller  values  of  nodal  displacement. 

In  order  to  obtain  an  improved  solution  for  the  elastica  problem,  the 
finite  element  mesh  was  revised  to  consist  of  10  B22  elements.  Due  to 
numerical  sensitivity  of  this  problem,  the  analysis  was  carried  out  in  two 
stages;  namely,  pre-buckling  and  post-buckling  stages.  For  the  pre-buckling 
stage,  the  load  control  option  with  direct  incrementing  procedure  was 
exercised.  The  load  was  increased  from  0  to  a  critical  value  Per  =  -  4.8  lb. 
For  the  post-buckling  stage,  the  analysis  was  switched  to  displacement  control 
with  an  automatic  incrementing  procedure.  In  both  cases,  a  tight  tolerance, 
i.e.  PTOL  =  MTOL  =  .05  was  specified. 
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Table  6.1.  Deformation  Modes  of  Elastica  vs.  Solution  Parameters. 


PTOLflbT'V^^ 

MTOL(ft-lbTV^ 

1 

5 

10 

20 

1 

10 

100 

Ck]  NUM;  User  suggested  number  of  increments 


MTOLT 
PTOL  J 


User  defined  moment  and  force  tolerances  for 
accuracy  control  (notice  that  the  max.  value  of  end 
force  for  this  problem  is  about  35  lb.) 


kt. 
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The  load-deflection  response  obtained  from  ABAQUS  is  compared  with  the 
closed  form  solution  [31]  as  shown  in  Fiy.  6.2  Obviously,  excellent 
agreement  was  found.  The  deflected  shapes  of  the  elastica  corresponding  to 
different  loads  were  plotted  in  Fig.  6.3.  These  shapes  appear  to  be  in 
reasonable  agreement  with  those  in  reference  [31]. 

6.2  An  Elastic-Plastic  Plate  Subjected  to  Biaxial  Loading. 

A  square  plate  subjected  to  biaxial  stresses  was  considered.  The 
primary  objective  of  this  benchmark  is  to  test  the  numerical  stability  of 
ABAQUS  under  elastic-plastic  loading,  unloading  and  neutral  loading.  The 
material  was  assumed  to  bi linearly  elastic-plastic  with  kinematic  strain 
hardening.  The  material  constants  are: 

Young's  modulus:  E  =  30  x  10^  Ksi 
Poisson's  ratio:  v  =  0.3 
Initial  yield  stress:  Oq  -  30  ksi 
Plastic  modulus:  Ep  =  1.5  x  10^  ksi 

The  loading  path  was  specified  in  3  successive  steps: 

Step  1.  0  <  Ox  1 

ay  =  0 

Step  2.  Ox  reduced  from  40.5  ksi  to  Oq  (unloading) 

Oy  =  0 

Step  3.  Ox  and  Oy  were  varied  according  to 
Ox^  +  oya  -  Ox  ®y  ®  ksi)^ 

Until  Oy  =  0  and  Ox  =  -  40  ksi 
The  above  loading  path  is  depicted  in  Figs.  6.4.  and  6.5. 
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Fig.  6.2.  Load  -  Deflection  Response  of  An  Elastica 
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Since  the  plate  was  subjected  to  a  homogeneous  state  of  stress,  only 
one  plane  stress  element,  CPS4,  was  used.  Initially,  direct  procedure  with  a 
loading  increment  aP  =  40000  psi  was  specified.  Although  ABAQUS  gave 
convergent  solution,  the  solution  deviates  from  the  analytical  result  [33] 
as  indicated  in  Fig.  6.6.  Analysis  was  also  done  by  using  ADINA.  Solution 
convergence  could  not  be  reached  as  soon  as  the  plasticity  was  initiated. 

In  order  to  find  a  suitable  loading  increment,  the  automatic  load 
incrementation  option  of  ABAQUS  was  exercised  for  7-load  increments.  The 
solution  correlates  with  the  analytical  result,  as  shown  in  Fig.  6.6,  very 
closely.  Once  a  convergent  load  increment  was  found  from  the  automatic 
incrementation  option,  a  direct  control  was  used  to  complete  the  analysis. 
Again,  as  seen  in  Figs.  6.6  and  6.7,  the  stress  and  strain  printout  from 
ABAQUS  are  indistinguishable  from  those  of  the  analytical  solution.  The 
analysis  was  also  conducted  by  using  ADINA  with  the  same  loading  increments. 
The  strain  output  from  ADINA  are  slightly  different  from  those  of  ABAQUS  as 
shown  in  Fig.  6.8. 

From  the  experience  of  this  example,  it  is  felt  that  the  automatic 
step  control  is  t..i  extremely  desirable  feature  for  searching  for  convergent 
load  increment  in  nonlinear  analysis. 

6.3  A  Cantilever  Beam  Subjected  to  an  End  Force 

Large  displacement  of  a  cantilever  beam  was  considered.  The  beam  was 
subjected  to  a  concentrated  force  at  the  free  end  (shown  in  Fig.  6.9)  and  two 
load  cases  were  examined:  Case  a)  Conservative  load  -  load  always  remained 
vertical,  and  ca  e  b)  Nonconservative  load  -  load  remained  normal  to  the 
deformed  neutral  axis  of  the  beam.  The  beam  was  assumed  to  have  a  rectangular 


Total  Strain  (  %) 
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Direct  Increment  with  Initial  AP  =  a 

0 

Direct  Increment  with  Initial  AP  =  40  ksi  >  o 


Fig,  6.6.  Strains  vs.  Solution  Increments  of  Biaxially  Loaded  Plate 


Fig.  6.7.  Plastic  Strains  vs.  Solution  Increments  of 
Biaxially  Loaded  Plate 


Total  Strain  (  %) 
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-  ABAQUS 

_  ADINA 


Fig.  6.9.  A  Cantilever  Bean  Subjected  to  a  Concentrated  Force. 


cross  section  with  linearly  elastic  material.  The  finite  element  model 
consisted  of  10  B21  (beam)  elements  with  11  nodes. 

When  using  ABAQUS,  the  automatic  load  control  with  fairly  large 
convergence  tolerances  was  specified,  i.e. 

PTOL  =  15%  •  (Max.  Reaction  Force) 

MTOL  =  2%  •  (Max.  Reaction  Moment) 

The  horizontal  and  vertical  displacements  at  the  free  end  for  the  two  load 
cases  considered  were  plottec  in  Figs.  6.10  -  6.12  against  the  solution 
obtained  by  Horrigmoe  [34], 

Due  to  the  difference  in  convergence  characteristics,  the 
non-conservative  load  case  required  much  greater  number  of  solution 
increments,  thereby  more  CPU  time,  than  the  conservative  load  case.  A 
comparison  is  given  as  the  following: 

No.  of  Increments  Total  Iteration  Cycles 


Conservative  Load 

10 

30 

Non-Conservative  Load 

30 

115 

It  was  somewhat  surprising  that  the  two  load  cases  gave  such  remarkable  - 
difference  in  convergence. 

The  same  problem  was  also  run  using  ADINA.  Since  ADINA  cannot  handle 
deformation-dependent  loading,  only  the  conservative  load  case  was  considered. 
Due  to  the  stiffening  load-deflection  response  of  this  problem,  solution 
convergence  was  extremely  difficult  to  achieve  when  the  MNR  method  with 
equilibrium  iterations  was  optioned.  After  several  unsuccessful  trials 
finally  a  procedure  with  800  load  increments  together  with  equilibrium 
iterations  was  necessary  to  obtain  a  convergent  solution  comparable  to 


Load  P  (KN) 
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Horrigmoe's  result.  On  the  other  hand,  an  ADINA  run  with  100  load  increments 
and  no  equilibrium  iterations  was  carried  out  and  the  result  is  similar  to 
that  of  800-increment  analysis.  This  phenomenon  illustrated  the 
unpredictability  of  the  MNR  algorithm  when  it  is  applied  to  a  nonlinear 
structural  problem  with  stiffening  effect. 

It  is  noted  that  both  ABAQUS  and  ADINA  results  agree  quite  well;  they 
are  however  slightly  different  from  those  obtained  by  Horrigmoe  and  Bergan 
[34],  especially  the  horizontal  component  of  the  displacement  at  the  free  end 
(Fig.  4.10). 

6.4  Large  Deformation  of  a  Spherical  Cap 

A  spherical  cap  subjected  to  a  concentrated  load  at  the  crown  as  shown 
in  Fig.  6.13  was  analyzed  using  both  ABAQUS  and  ADINA.  Again,  the  purpose  of 
this  example  is  to  demonstrate  the  relative  merits  between  the  full  and 
modified  Newton -Raphson  algorithms. 

The  material  of  spherical  cap  was  assumed  to  be  linearly  elastic.  This 
problem  was  previously  analyzed  by  Stricklin  [35],  Mescal  1  [36],  and  Bathe  in 
ADINA  manual  [11].  As  seen  in  the  load  -  deflection  plot  (Fig.  6.13), 
initially  the  cap  exhibits  softening  behavior.  After  the  apex  is  deformed  to 
a  point  passing  the  baseline  through  the  two  supports,  the  bending  action  in 
the  cap  is  transformed  into  the  membrane  action,  and  thus  stiffening  of  the 
structure  can  be  observed. 

The  cap  was  modeled  by  ten  8-node  axi symmetric  elements,  i.e.  CAX8R 
elements.  To  run  ABAQUS,  an  automatic  procedure  with  the  following  control 
parameters  was  used: 
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P 


R  .  4.76 
h  =  0.01576  in 
H  =  0.0859  in 
0  s  10.9“ 

\  *  6 

E  »  10  X  10*  Ib/in* 

V  z  0.3 

p  =  0.000245  lb“sec*/in 


Fig.  6.13.  Load-Deflection  Curves  for  Spherical  Shell 
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Suygested  Number  of  Increments:  NUM  =  5 

Force  Tolerance:  PTOL  =  100  {  “  1%  of  Max.  Reaction) 

The  analysis  from  ABAQUS  was  completed  with  a  total  of  36  solution  increments 
(including  iterations),  and  the  results  are  shown  in  Fig.  6.13. 

When  ADINA  was  used  to  run  the  same  problem,  considerable  difficulty  in 
solution  convergence  was  experienced  if  equilibrium  iterations  were  exercised. 
Even  with  different  specifieii  values  of  convergence  tolerance  and  maximum 
iteration  cycles,  the  attalyses  were  interrupted  for  loads  between  20  -  50  lb. 
After  several  attempts,  a  procedure  with  80  load  increments  and  no  iterations 
was  successful  to  obtain  a  solution  agreeable  to  the  known  results.  Once 
again,  this  example  demonstrated  the  numerical  difficulty  of  the  MNR  algorithm 
whereas  the  FNR  method  appears  to  be  superior. 

6.5  Buckling  of  A  Spherical  Shell 

In  general,  ABAOUS  cannot  be  used  for  post-buckling  analysis,  since  the 
program  does  not  have  a  special  solution  algorithm  for  handling  such  case. 
Nevertheless,  post-buckling  analysis  can  be  performed  for  structures  subjected 
to  simple  loadings,  such  as  a  concentrated  force  applied  at  the  point  of 
symmetry.  For  this  purpose,  a  displacement  control,  as  opposed  to  the  load 
control,  must  be  specified  for  the  intended  analysis. 

A  spherical  shell  with  linearly  elastic  material  properties  was 
selected  for  buckling  anaysis.  The  purpose  of  such  analysis  is  to  determine 
the  convergence  characteristics  of  shell  elements  in  ABAQUS  when  significant 
change  of  geometry  is  involved.  The  same  problem  has  previously  been  analyzed 
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by  several  researchers  [24,  34,  37]  due  to  its  simplicity  in  geometric  shape 
and  dramatic  post-buckling  features.  Dimensions  of  material  constants  of  the 
shell  are  given  as  below; 

Radius  R  =  2450  mm 

Side  Dimension  a  =  754.9  nm 

Thickness  of  Shell  h  =  99.45  mm 

Young's  Modulus  E  =  68.95  N/mm^ 

Poisson's  Ratio  v  =  0.3 

The  shell  is  supported  by  hinges  in  such  a  way  that  it  can  rotate  freely  along 
the  tangent  direction  to  the  curved  edges,  but  no  translation  is  allowed. 

With  this  boundary  condition,  coordinate  transformation  must  be  introduced  in 
using  ABAQUS,  since  the  constrain  conditions  are  specified  in  reference  to  the 
local  coordinate  at  the  boundary  nodes. 

Taking  advantage  of  symmetry,  only  one-quarter  of  the  shell  surface  was 
modeled  by  S8K  shell  elements.  Three  different  meshes  were  considered  for  the 
intended  analysis;  namely  2  x  2  ,  3  x  3  and  4x4  meshes.  The  load  vs. 
maximum  deflection  responses  of  the  three  meshes  are  shown  in  Fig.  6.14.  It 
is  seen  that  the  load-deflection  responses  obtained  from  the  3x3  and  4x4 
meshes  are  identical;  whereas  the  prebuckling  response  of  2  x  2  mesh  is 
similar  to  those  of  finer  meshes,  but  the  post-buckling  response  is  far 
different.  In  fact,  when  the  2x2  mesh  was  loaded  near  the  critical  value, 
the  shell  suddenly  jumped  to  a  new  equilibrium  state  (membrane  action), 
bypassed  the  lower  post-buckling  range.  However,  when  the  load  was  reduced 
from  the  membrane  state,  the  lower  post-buckling  curve  was  traced,  but  quite 
different  from  those  of  the  finer  meshes.  The  implication  of  this  phenomenon 


3x3  and  4x4  meshes 
2x2  mesh 


Response  Due  to  Force  Control 


100  150  200  ^ 

Deflection  (mm) 


Post-Buckling  Behavior  of  a  Spherical  Shell 
With  Different  Mesh  Sizes 
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is  that  a  finite  element  mesh  is  sufficient  for  pre-buckling  analysis,  but 
may  not  be  the  case  for  representing  post-buckling  behavior. 

For  comparison,  the  load  deflection  response,  obtained  from  ABAQUS 
using  4x4  mesh,  was  plotted  in  Fig.  6.1b  togehter  with  those  of  Leicester 
[37]  and  Horrigmoe  and  Beryan  [34].  Small,  but  significant  deviation  can  be 
seen  between  ABAQUS  results  and  others.  Such  deviation  might  have  been  caused 
by  the  use  of  TRANSFORM  option  in  ABAQUS  for  defining  the  boundary  restraints 
with  respect  to  the  local  coordinates  along  the  curved  edge  of  the  shell. 

6.6  A  Centrally  Cracked  Plate 

A  rectangular  plate  containing  a  plane  crack  at  its  center  was 
considered  as  another  example  (Fig.  6.16).  The  plate  is  subjected  to  uniform 
pressure  Py  at  its  two  ends.  The  material  properties  were  assumed  to  be 
elastic  perfectly  plastic  ob.jying  von  Mises  yield  criterion  and  the  material 
constants  are 

Young's  Modulus:  E  =  27.4  x  10^  ksi 

Poisson's  Ratio:  v  =  0.3 

Yield  Stress:  oq  =  10  ksi 

The  purposes  of  this  benchmark  are  two  fold:  1)  to  check  the  correctness  of 
3/0  solid  element  in  ABAQUS  for  elastic-plastic  analysis,  and  2)  to  determine 
the  sensitivity  of  solutions  corresponding  to  different  specified  force 
tolerance  values.  The  same  problem  was  previously  anayzed  in  [38]. 

The  plate  was  represented  by  both  2/0  (CPE8)  and  3/0  (C3020)  solid 
elements,  respectively.  From  symmetry,  only  one-quarter  of  the  plate  was 
modeled.  For  the  3/0  elements,  three  layers  of  nodes  were  used.  To  utilize 
the  node  and  element  generation  features  in  ABAQUS,  the  node  number  starts  at 
the  crack  tip  in  the  plane  of  symmetry  and  increases  circumferentially 


Deflection  (nin) 


6.15.  Load  vs.  Deflection  at  the  Center  of  Shell 
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Fig,  6.16.  A  Centrally  Cracked  Plate 
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outward  (Fig.  6.17).  All  nodes  are  restrained  from  movement  in  the  thickness 
direction  so  that  a  plane  strain  state  prevails.  The  same  mesh  and  similar 
boundary  conditions  were  used  for  the  2/D  model. 

In  running  this  problem,  the  following  solution  parameters  were 
specified  for  both  2/D  and  3/D  analyses: 

a)  Automatic  loading  increment 

b)  Maximum  number  of  increments:  INC  =  50 

c)  Maximum  number  of  Iteration  cycles  in  an  increment: 

CYCLE  =  6  (Default  value) 

d)  Suggested  number  of  increments:  NUM  =  5 

e)  Force  tolerance  :  PTOL  =  varied. 

Several  runs  were  made  for  both  2/D  and  3/D  models  by  varying  the  force 
tolerances  to  control  solution  convergence.  The  values  of  tolerance  (PTOL) 
used  are  specified  in  Table  6.2. 

Some  of  the  analysis  results  corresponding  to  the  above  runs  were 
plotted  against  the  load  for  comparison,  i.e. 

a)  Load  vs.  vertical  deflection  of  node  #41  (center  of  crack)  - 
Figs.  6.18  and  6.19. 

b)  Vertical  normal  stress  oy  near  crack  tip  (xg  =  0.1“  ,  =  0.07", 

from  the  tip)  -  Figs.  6.20  and  6.21. 

c)  Horizontal  stress  (thickness  direction)  at  the  same  location  as 
b  -  Figs.  6.22  and  6.23. 

It  is  seen  from  Fig.  6.18  that  the  vertical  component  of  displacement  at  node 
#41  obtained  from  the  2/0  model  was  not  sensitive  to  the  convergence 
tolerance.  The  same  phenomenon  can  be  said  for  the  stresses  (Figs.  6.20  and 
6.22).  On  the  other  hand,  the  results  obtained  from  the  3/0  model  for  the 


I 


12.6 


Vertical  Displacerrent  (in.) 

0.16  0.214  0.32  0.40  0.46 


10.00  20.00 
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Table  6.2  Different  Tolerance  Values 
for  the  Crack  Problem 


Run  No. 

Description 

Force 

Tolerance 

1 

2/D  Model 

350 

lb  (.42%) 

2 

2/D  Model 

3,500 

(4.25%) 

3 

3/D  Model 

3,500 

(0.17%) 

4 

3/D  Model 

8,000 

(0.4%) 

5 

3/0  Model 

15.000 

(0.75%) 

6 

3/0  Model 

35,000 

(1.7%) 

♦Note; 

Percentage  for  force  tolerance  was  calculated  on  the  basis  of  maximum 
reaction  force  (absolute  value)  at  supports. 
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same  problem  is  rather  sensitive  to  the  tolerance  values  (note  that  the  value 
PTOL  =  35,000  lb  is  only  1.7%  of  the  maximum  reaction  force,  well  within  the 
recommended  range  stated  in  ABAQUS  manual).  As  the  plastic  zone  was  further 
spread  out  around  the  crack  tip,  the  errors  for  near  crack  tip  stresses  are 
becoming  even  greater  and  erratic. 

Another  point  o’  consideration  is  the  computation  (CPU)  time  required 
for  solving  a  3/0  elastic-plastic  crack  problem.  The  CPU  time  increases 
significantly  with  the  user's  assigned  convergence  tolerance.  A  comparison  of 
CPU's  for  different  cases  of  3/0  runs  is  given  in  Table  6.3.  Since  the  runs, 
except  the  case  5,  were  not  completed  at  100%  of  the  intended  loading,  the 
CPU  time  for  completing  each  problem  was  estimated  in  the  last  column  of  Table 
6.3.  It  is  seen  clearly  that  the  computer  time  required  for  analyzing  this 
relatively  small  nonlinear  3/0  problem  is  in  the  order  of  1.5  -  4  hrs  IBM  CPU! 

The  iteration  h'stories  and  automatic  load  cut-back  in  ABAQUS  for  the 
3/0  analyses  are  shown  in  Fig.  6.24.  It  is  seen  that  more  frequent  load 
cut-back  occurs  for  smaller  force  tolerances,  and  thereby  it  requires  more 
iteration  cycles  to  complete  the  problem. 

From  the  experience  of  this  example,  the  following  comments  are  drawn: 

i)  Analysis  results  for  3/0  problems  are  sensitive  to  the 

convergence  tolerance  specified  by  the  user.  This  phenomenon  is 
especially  acute  for  stresses  (or  strains)  as  opposed  to  the 
displacements.  It  is  noted,  however,  that  this  numerical  difficulty 
did  not  appear  for  the  2/0  model. 

ii)  For  application  problems,  it  is  rather  difficult  for  the  user  to  come 
up  with  a  tolerance  value  (PTOL  or  MTOL)  which  would  arrive  at  a 
reasonable  solution  accuracy  and  computation  time. 
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Table  6.3  CPU  Time 


vs.  Tolerance  of  3/D 

Crack  Model 

Case 

PTOL 

CPU 

Spent 

No.  of 
Iterations 
Compl eted 

Ave.  CPU 
Per 

Iteration 

%  of  Total 
Load 

Completed 

Estimated 

CPU 

3 

3,500  lb 

120  min 

.  26 

5.30  min 

67% 

230  min. 

4 

8,000 

180 

39 

4.64 

87.9% 

205 

5 

15,000 

130 

28 

4.62 

100% 

130 

6 

35,000 

90 

17 

4.62 

88.7% 

95.3 

Load  (%) 
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®-«-0  PTOL  =  3500  lb 
^  PTOL  =  8000  lb 

*  *-  »  PTOL  =  15000  lb 
X-X-X  PTOL  =  35000  lb 


o 

rsi 


Fig.  6.24.  Iteration  History  vs.  Automatic  Load  Cut-Back 
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6.7  Dynamic  Response  of  A  Spherical  Cap 

A  spherical  cap  subjected  to  uniform  step  pressure  was  considered  as  a 
benchnmark  to  test  the  dynamic  algorithms  being  used  in  ABAQUS.  The  same 

problem  was  previously  analyzed  by  Bathe,  Ramm  and  Wilson  [39].  For  the 

present  study,  the  following  three  cases  were  considered: 

Case  1.  Linearly  el.istic  material  and  small  deformation 
Case  2.  Elastic-plastic  material  and  small  deformation 

Case  3.  Elastic-plastic  material  and  large  rotations 

The  primary  purpose  of  this  benchmark  is  to  study  the  sensitivity  of  solution 
convergence  to  the  parameters  as  below: 

i)  Size  of  time  increment  At 

ii)  Artificial  dampinj  coefficient  a 

iii)  Convergence  tolerances  PTUL  and  HAFTOL. 

In  addition  to  ABAQUS  runs,  analyses  were  also  performed  using  ADINA  for  the 
purpose  of  comparing  results. 

For  case  1,  i.e.  linearly  elastic,  small  deformation  analysis,  the 
material  constants  are 

Young's  modulus  F.  =  10.5  x  10^  psi 

Poisson's  ratio  u  =  0.3 

Density  p  =  2.45  x  lO"**  lb.  sec^/in 

The  solution  parameters  used  for  ABAQUS  and  ADINA  are  shown  in  Table  6.5  by 
10-axi symmetric  elements  and  63  nodes.  The  results  obtained  from  both  ABAQUS 
and  ADINA  are  plotted  in  Fig.  6.25.  Solutions  from  the  two  programs  agree 
quite  closely.  Although  the  total  number  of  time  increments  required  for 
ABAQUS  is  only  one-half  of  that  needed  by  ADINA  in  order  to  achieve  the  same 


A 
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Table  6.4  Solution 

ABAQUS 

Time  integration  Hilber-Hughes  operator 

a  =  -.05 

Convergence  control  Residual  Force: 

HTOL  =  4000  lb. 

Time  Increment  At  =  10”4  sec. 


Parameters 

ADINA 

Newmark  operator 
a  a  0 

Residual  Energy 
ETOL  •  10-3 

At  a  0,5  X  10*3  sec. 


Deflection  at  Apex  Wo  (in) 
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Fig.  6.25.  Comparison  of  Linear  Response  of  a  Spherical 
Cap  Between  ABAQUS  and  ADINA 
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order  of  numerical  accuracy,  the  actual  CPU  time  of  ABAQUS  for  this  small 
problem  is  more  than  five  times  that  of  AOINA.  Although  ABAQUS  appears  to  be 
slower  than  AOINA  in  handling  small  problems,  it  would  be  unfair  to  make  a 
definite  statement  in  regard  to  their  relative  numerical  expediency  for 
dynamic  analysis  due  to  the  difference  in  their  programming  structure  and 
intended  purpose  between  these  two  codes. 

By  varying  the  size  of  time  increment  At  but  leaving  other  parameters 
in  Table  6.6  unchanged,  additional  linear  dynamic  analyses  were  conducted  to 
determine  the  solution  sensitivity  of  ABAQUS  and  AOINA  with  respect  to  At. 

The  deflection  histories  at  the  apex  of  the  cap  were  plotted  in  Fig.  6.26  for 
ABAQUS  and  Fig.  6.27  for  AOINA,  respectively.  It  is  seen  that  ABAQUS  solution 
starts  to  deteriorate  when  At  was  increased  to  5  x  10"^  sec;  whereas  the 
solution  from  AOINA  starts  to  deviate  for  At  =  2  x  10"5  sec.  This  comparison 
indicates  that  the  numerical  algorithm  for  dynamic  analysis  in  AOINA  is 
slightly  more  sensitive  to  the  increment  size  At  than  that  in  ABAQUS. 

The  second  parameter  to  be  investigated  is  the  use  of  artificial 
damping  coefficient  o  in  ABAQUS.  Three  different  damping  values  were  chosen, 
namely  a  =  0,  -.05,  -1/3,  and  the  linear  analyses  were  made  with  a  constant 
time  increment  At  =  10"^  sec  and  the  direct  incrementing  procedure.  The 
displacements  at  the  apex  vs.  time  are  shown  in  Fig.  6.28  and  the  damping 
coefficient  did  not  seem  to  cause  any  significant  difference  in  results.  This 
is,  of  course,  expected  since  At/T  (T  =  fundamental  period)  ratio  is  less  than 
1/40,  well  below  the  range  where  numerical  dissipation  may  occur  [12], 

However,  the  damping  ratio  appears  to  have  a  profound  effect  on  automatic 
incrementing  procedure  in  ABAQUS.  For  example,  when  the  same  input  data  were 
run  by  switching  from  direct  to  automatic  procedure,  solution  convergence 
becomes  extremely  difficult  for  a  =  0.  For  this  case,  ABAQUS  cut  back  the 

i 

I 


Deflection  At  Apex  Wo  (in) 
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Fig.  6.27.  Linear  Dynamic  Response  of  a  Cap  Using  ADINA 
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suggested  time  increment  At  =  lO'^  sec  to  10  -  subincrements  so  that  the 
calculated  half-step  residue  force,  R(l/2),  was  within  the  assigned  tolerance. 
The  changes  in  R(l/2)  and  the  vertical  deflections  at  the  apex  of  the  cap  vs. 
time  during  the  first  increment  cut-back  are  shown  in  Table  6.5.  One  can  see 
that  the  maximum  R(l/2)  at  the  initial  step  was  2.56  x  10®  lbs,  much  greater 
than  the  assigned  tolerance  20,000  lbs.  Figs.  6.29-31  show  the  plots  of  the 
maximum  R(l/2)  vs.  time  for  a  =  0,  -0.05,  -1/3,  respectively.  It  is  clear 

that  the  calculated  R(l/2)  is  very  sensitive  to  the  damping  a.  On  the  one 
hand,  if  a  zero  or  small  damping  ratio  is  used,  the  number  of  time  increments 
required  for  automatic  procedure  in  ABAQUS  would  be  enormous.  On  the  other 
hand,  if  a  large  damping  ratio  is  used,  ABAQUS  will  increase  the  size  of  At 
automatically,  thereby  introducing  significant  numerical  errors  for  the  long 
time  response  as  seen  in  Fig.  6.32,  for  a  =  -1/3.  The  manner  of  step  changes 
(i.e.  time  vs.  increment  numbers)  is  also  shown  in  Fig.  6.33.  Therefore,  as 
suggested  in  ABAQUS  manual,  the  best  damping  ratio  appears  to  be  a  *  -0.05 
whf-n  the  automatic  incrementing  is  exercised.  However,  for  direct  procedure 
the  value  of  damping  ratio  has  no  si<inificant  effect  on  the  analysis  results 
so  long  as  At/T  <  .1, 

For  the  case  2,  eliistic-plastic  material  with  small  deformation  was 
considered  and  the  material  constants  are 
Initial  yield  stress  oq  =  24  ksi 
Plastic  modulus  Ep  =  0.21  x  10^  ksi 

(Isotropic  strain  hardening  rule) 

Solution  parameters  used  for  ABAQUS  are 
aT  =  10"®  sec. 

o  =  -0.05 


R(l/2)  =  20,000  lbs. 
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Table  6.5  Change  of  HAFTOL  and  Vertical  Displacement 

at  the  Apex  During  the  First  Increment  Cut-Back 

T(*)  DT  R{l/2)  U 


l.OOOE-5 

1.56E-7 

2560000 

(**) 

1.560E-6 

8.8bE-7 

706 

-1.569E-6 

1.040E-6 

II 

27100 

-1.266E-5 

2,026E-6 

11 

11500 

-2.129E-5 

2.822E-6 

II 

11500 

-3.055E-5 

3.917E-6 

l.llE-6 

31800 

-5.142E-5 

5.023E-6 

II 

27500 

-8.473E-5 

6.129E-6 

M 

8590 

-1.261E-4 

7.335E-6 

12500 

-1.653E-4 

8.616E-6 

1.38E-6 

33100 

-2.299E-4 

l.OOOE-5 

M 

33800 

-3.144E-4 

{*) 

T . 

Real  Time 

T(*) . 

Normalized  Time 

DT  . 

Time  Increment 

R(l/2)  .. 

Half  Residual 

(**) 

Cut  back 

due  to  excessive  residual 

Residue  (10 


Displacement  at  Apex  wo  (in) 
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Fig.  6.32.  Effect  of  Algorithmic  Damping  vs.  Stepping  Method 


Fraction  of  Total  Solution  Time 
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and  direct  procedure  was  exercised.  The  displacement  history  at  the  apex 
obtained  from  ABACUS  is  compared  with  that  of  ADINA  in  Fig.  6.34.  Both 
results  appear  to  agree  quite  closely. 

The  third  case  is  to  consider  both  the  geometric  (large  rotations)  and 
material  nonlinearities.  For  this  case,  both  the  automatic  and  direct 
incrementing  procedures  were  employed  and  results  are  plotted  in  Figs.  6.35  . 
Again,  ABAQUS  results  are  in  general  agreement  with  ADINA,  although  some  minor 
difference  can  be  seen.  It  is  noted,  however,  that  when  the  automatic 
procedure  was  used  for  ABAQUS,  some  unnecessary  increment  cut -back  was 
experienced  in  the  same  way  as  the  linear  case.  This  is  considered  to  be  the 
shortcoming  of  the  automatic  procedure.  Nevertheless,  the  reviewer  still 
believes  that  this  option  is  very  helpful  for  handling  a  nonlinear  problem  in 
which  case  the  convergent  time  increment  is  not  known  a  priori. 


ection  at 


Displacement  at  Apex  Wo  (in) 


C.8  Dynamic  Analysis  of  A  Semi -Spherical  Shell 


The  linear  dynamic  response  of  a  semi-spherical  shell  subjected  to 
uniform  pressure  was  considered  as  another  dynamic  case.  The  unique  part  of 
this  problem  is  that  its  dynamic  response  involves  a  wide  spread  of 
frequencies  and  mode  shapes.  The  motives  of  this  benchmark  are  two-fold: 
i)  To  compare  the  results  obtained  from  both  ABAljUS  and  ADINA  with  existing 
analytical  solution  [22],  and  ii)  to  examine  the  effect  of  numerical 
dissiption  when  a  range  of  vibration  modes  are  involved. 

The  benchmark  under  consideration  was  previusly  solved  by  Kraus  and 
Kalnins  [40].  The  geometry  and  coordinate  system  of  the  shell  are  shown  in 
Fig.  6.36.  The  general  solution  for  displacement  of  the  shell  was  given  in 
[39]  and  it  is  outlined  in  the  Appendix  A.  For  a  specific  geometry  and 
material,  i.e. 

h=lin  a=20in 

E  =  30  x  10^  psi  V  =  0.3  p  =  1  Ib-sec^/in** 

and  simply  supported  condition,  the  vertical  displacement  at  the  crown  has  the 
expression 


w 


1.33  X  10-^  Pg 


L 

i=l 


mi 


(1  -  cos  273.86  t] 


(1) 


where  Pg  is  the  applied  uniform  pressure;  m-j ,  the  modal  participation  factor; 

,  non-dimensional i zed  frequency  defined  in  Appendix  A.  In  order  to  obtain  a 
sufficiently  accurate  solution,  a  total  of  15-mode  shapes  are  required.  This 
problem  was  anayzed  by  both  ABAQUS  and  ADINA.  The  finite  element  model  used 
is  a  uniform  mesh  consisting  of  20  -  axisymmetric  elements  and  123  nodes. 
Finite  element  results  in  the  form  of  displacement  and  velocity  histories  at 
the  apex  are  compared  with  the  analytical  solution  as  seen  in  Figs.  6.37  -  38. 


Displacement  at  Apex  (10 


Velocity  at  Apex  (10’  In/sec) 
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It  is  seen  that  other  than  spurious  oscillations  and  a  small  amount  of  period 
elongation  from  ADINA  solution  (due  to  Wilson  e  operator),  all  solutions  are 
generally  in  good  agreement. 

Additional  runs  were  made  with  ABAQUS  by  varying  the  increment  size  At 
and  damping  ratio  a,  i.e., 

At  =  3  X  10"‘*  sec.  ,  6  X  10"^  sec.  ,  12  x  lO"**  sec. 

a  =  0.  ,  -  1/3 

As  seen  in  Figs.  6.39  and  6.40,  the  solution  is  slightly  affected  by  the 
increment  size,  but  rather  insensitive  to  the  damping  ratio  as  found  out  from 
the  previously  example. 

6.9  Longitudinal  Pulse  in  a  Prismatic  Bar 

For  dynamic  analysis,  ABAQUS  adopted  a  numerical  operator  called  "a  - 
method"  proposed  by  Hilber  and  Hughes  [12,  13],  This  operator,  as  described 
previously,  possesses  a  single  parameter,  namely  a,  which  is  adjustable  for 
controlling  numerical  damping.  With  an  appropriate  a-value,  the  operator 
damps  out  the  higher  frequency  components  and  preserve  all  the  necessary  low 
frequency  modes. 

To  test  the  numerical  performance  of  ABAQUS  in  dynamic  analysis, 
longitudinal  wave  propagations  in  a  prismatic  bar  are  considered.  Two 
different  cases  wore  analyzed  by  ABAQUS: 

Case  1  -  Elastic  wave  of  a  triangular  pulse  (Fig,  6.41a) 

Case  2  -  Elastic-plastic  wave  of  a  ramp  pulse  (Fig.  6.41b) 


Displacement  at  Apex  (10*  in) 
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Fig.  6.39.  Displacement  Response  vs.  Different  Time  Increments 
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Material  properties  of  the  bar  are 

E  =  10**  Ksi  .  Ep  =  2  X  103  Ksi 

oq  =  20  Ksi  ,  P  =  2.5  X  ID"**  Ib-sec^/in** 

To  prepare  the  finite  element  model,  some  consideration  must  be  given  in  the 

determination  of  mesh  size.  Basically,  there  should  be  sufficient  elements  to 
capture  the  profile  of  the  wave  front,  namely  the  rising  pulse.  Using  the 
elastic  wave  as  a  basis  of  calculation,  we  find 

Longitudinal  wave  velocity:  Cg  =  •/  E/p  =  w  x  lO^  in/sec. 

Rising  time  of  the  pulse:  tp  =  5  x  sec. 

Length  of  the  wave  front:  Lf  =  tp  Cg  =  1  in. 

We  used  2-elements  to  cover  the  wave  front,  hence  the  mesh  size  in  the  region 
of  propagating  wave  is  Ax  =  0.5  in. 

The  next  point  of  consideration  is  the  use  of  time  increment  At.  For 
wave  propagation  analysis.  At  was  suggested  to  be  in  the  range:  Tp/lO  £  At  £ 
Tn/2,  Tp  =  highest  period  of  the  system.  For  the  triangular  pulse,  we  found 
^n  "  ^200  =  5  X  10"®  sec.  Therefore,  a  At  =  .5  x  10“®  sec.  was  selected  for 
the  analysis. 

After  the  mesh  size  and  lime  Increment  were  determined,  solution  of 
this  simple  wave  problem  from  AI  AQUS  was  rather  straightforward.  Fig.  6.42 
shows  the  elastic  stress  pulse  in  the  bar  for  t  =  SOps.  Elastic-plastic  wave 
propagations  are  shown  in  Figs.  6.43-46  for  instants  t  =  30,  60,  90  and  120  p 
sec,  respectively.  It  is  seen  that  A8AQUS  gave  almost  identical  results  as 
the  analytical  solutions  (Appendix  B). 

As  we  consider  the  results  of  Case  2,  the  elastic-plastic  wave  consists 
of  three  parts:  elastic  wave  front  (Fig.  6.43)  travelling  at  a  speed  Cg, 
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Fig.  6.42.  Elastic  Stress  Pulse 


plastic  wave  front  at  speed  Cp  (Fig.  6.44)  and  elastic  wave  back  due  to 
plastic  unloading  (Fig.  6.45).  At  a  certain  point,  the  wave  back  begins  to 
cut-off  the  plastic  front  and  causes  a  secondary  wave,  as  seen  in  Fig.  6.46, 
which  induces  a  stress 

S  =  ■t('^  -l)(a  -a) 

where  o  =  50  ksi;  a,  stress  at  an  instant  t  after  plastic  unloading.  The 
secondary  wave  propagated  both  forward  and  backward  at  a  speed  Cg. 

The  same  problem  was  previously  analyzed  by  Hartzman  and  Hutchinson 
[23]  by  using  2/0  4-node  elenents  and  explicit  integration  scheme.  In  their 
solutions,  considerable  spur'ous  oscillations  were  indicated  in  both  elastic 
and  elastic-plastic  waves.  From  ABAQUS  runs,  no  spurious  inodes  were  seen  when 
a  small  damping  ration,  a  =  -0.05,  was  imposed.  The  CPU  time  required 
for  solving  this  small  nonlinear  dynamic  problem  was  quite  excessive  on  the 
IBM  370-158  : 

Case  I  =  25  min.,  Case  2  =  17  min. 

6.10  Execution  Time 

One  question  that  wou  d  have  been  asked  by  many  users  is  how  efficient 
is  ABAQUS,  say,  as  compared  :o  AOINA?  It  is  difficult  to  answer  this  question 
without  some  sort  of  qualified  statement. 

First  of  all,  one  mus :  keep  in  mind  that  the  computational  efficiency 
of  a  finite  element  software  is  affected  by  several  factors.  These  include: 

i)  Program  architecture  and  coding  style 

ii)  Numerical  algorithm  adopted 

iii)  Application  problem  size 

iv)  Hardware  environment 
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Elastic-Plastic  Pulse,  t  =  60  u  Sec. 
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Elastic-Plastic  Pulse,  t  =  90 
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Fig.  6.46 
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First,  since  ABAQUS  was  developed  using  a  data  base  management  procedure  for 
the  purpose  of  coding  reliability,  invariably  some  degree  of  overhead  must 
be  paid  whenever  running  a  problem.  This  overhead  may  become  insignificant 
for  large  appliction  problems,  but  it  is  a  significant  portion  of  the  total 
CPU  time  for  small  size  problems.  Secondly,  the  nonlinear  solution  algorithm 
adopted  in  ABAQUS  is  b^sed  on  the  full  Newton-Raphson  method.  Although  this 
method  gives  better  convergence  characteristic  as  compared  to  the  modified 
Newton-Raphson,  this  advantage  is  often  off-set  by  its  excessive  computing 
effort.  Based  on  the  above  two  points,  it  appears  that  ABAQUS  would  be  less 
computationally  efficient  than  ADINA.  This  observation  was  also  agreed  upon 
by  ABAQUS  developers.  However,  it  was  argued  that  ABAQUS  becomes  more  CPU 
competitive  for  large  size  problems.  The  question  still  remains  that  how 
large  a  problem  has  to  be  in  order  for  ABAQUS  to  gain  its  CPU-competitiveness? 

To  answer  the  above  question,  two  benchmark  problems  were  run  using 
both  ABAQUS  and  ADINA:  i)  Elastic-plastic  analysis  of  a  cantilever  beam  and 
ii)  Elastic  response  of  a  Rubik  Cube.  Each  of  the  two  problems  is  described 
as  below. 

i)  A  Cantilever  Beam 

The  cantilever  beam  was  modeled  by  two  rows  of  8-node  plane  stress 
elements  as  shown  in  Fig.  6.47  The  number  of  elements  were  increased 
progressively  along  the  axis  of  the  beam,  thereby  the  number  of  nodes  were 
increased  accordingly.  Since  only  two  rows  of  elements  are  involved  for  this 
problem,  the  wave  front  for  ABAQUS,  or  bandwidth  for  ADINA  is  rather  small, 
even  though  the  total  number  of  degrees  of  freedom  can  be  expanded. 

The  beam  was  subjected  to  a  concentrated  force  at  its  free  end  and  the 
material  constants  used  for  the  analysis  are 


Young’s  modulus 
Poisson's  ratio 


E  =  30  X  103  Ksi 


V  =0.3 

Initial  yield  stress  oq  =  30  Ksi 
Plastic  modulus  Ep  =  1500  Ksi 

Both  the  material  and  geometric  nonlinearities  were  optioned  with 
20-load  increments.  Plasticity  was  initiated  during  the  second  load 
increment.  For  ABAQUS,  direct  incrementing  procedure  was  used.  The  total 
number  of  elements  varied  from  100  to  400.  Correspondingly,  the  number  of 
degrees  of  freedom  varied  between  810  to  3210.  For  this  problem,  the  CPU  time 
on  IBM  370-158  computer  for  both  ABAQUS  and  ADINA  is  plotted  against  the 
number  of  elements  as  shown  in  Fig.  6.48.  The  CPU  for  ABAQUS  increases 
linearly  with  the  model  size.  For  ADINA,  two  different  sizes  of  workspace 
(master  A-array)  were  specified;  namely  45  K  and  90  K  (words).  It  is  seen 
that  the  efficiency  of  AOINA  strongly  depends  on  the  size  of  workspace 
allocated.  This  is  due  to  the  fact  that  the  code  uses  an  out-of-core  blocking 
technique  in  which  the  structural  stiffness  is  partitioned  into  several  blocks 
and  a  great  deal  of  CPU  is  consumed  in  transferring  data  from  the  fast  core  to 
low-speed  disk  or  tape  units  or  vice  versa.  Obviously,  ADINA  is  quite 
efficient  for  small  size  models.  As  the  number  of  elements  was  increased, 
ADINA  loses  its  advantage  over  ABAQUS.  One  way  to  improve  ADINA' s  efficiency 
is  to  increase  the  core  size  for  execution.  Even  with  a  virtual  memory 
system,  there  is  always  an  upper  limit  on  the  maximum  core  size  that  can  be 
allocated. 

i  i )  A  Rubi k  Cube 

Differing  from  the  preceding  benchmark,  both  the  number  of  degrees  of 
freedom  and  its  bandwidth  increase  as  the  mesh  size  is  expanded.  The  cube  is 
subjected  to  a  concentrated  force  at  the  center  and  symmetric  boundary 
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Fig.  6.48.  Comparison  of  CPU  Time  Between  ABAQUS  and  ADINA  for  a  2/D  Problem 


Table  6.6  ADINA  Runs  for  Elastic-Plastic  Beam  Problem 


Problem 

Size 

Workspace 

Number 
of  Blocks 

CPU  for 

Stiffness 

Formation 

CPU  for 
Triangu- 
larize 

Step-by- 

step 

Total 

Total 

CPU 

50 

A5,000 

1 

0* 

0* 

4min  54  sec 

5fnin  3i,95sec 

50 

90,000 

1 

1  sec 

0 

4min  42  sec 

5"'if’  20,69sec 

100 

45,000 

3 

1  sec 

0 

limin  56sec 

12"'i'i35.99sec 

100 

90,000 

1 

0 

0 

7  min  42sec 

8"’ii  3.36sec 

150 

45,000 

8 

1  sec 

0 

25min  j  sec 

25'"in41 .28sec 

150 

90,000 

1 

0 

0 

limin  49Sec 

12mini2.39sec 

200 

90,000 

3 

1  sec 

0 

22niin  48Sec 

23f"if’30.5isec 

250 

90,000 

5 

1  sec 

0 

39min  45sec 

40'"if’10.8isec 

300 

90,000 

8 

1  sec 

0 

50min  54Sec 

5T"i"21.19sec 

350 

90,000 

19 

2  sec 

0 

Bomin  i5sec 

0.23sec 

350 

1,000,000 

1 

1  sec 

0 

34min  24sec 

35"'if’  8.73sec 

*Less  than  .1  second 


conditions  were  enforced.  Since  full  nonlinear  analysis  of  this  problem  is 
prohibitively  expensive,  only  linearly  elastic  response  was  considered.  The 
mesh  size  was  varied  from  two  to  five  divisions  in  all  three  directions  and  a 
mesh  of  5  X  5  X  5  is  shown  in  Fig.  6.49.  The  workspaces  specified  for  ABAQliS 
were  50  K  bytes  and  300  K;  for  ADINA,  90  K  and  300  K.  The  CPU  time  vs.  mesh 

size  is  plotted  in  Fig.  4.  Again,  for  small  problems  AOINA  is  more  efficient 

than  ABAQUS.  However,  as  the  size  of  the  problem  was  increased  to  be  5  x  5  x  5 
i.e.  125  elements  and  56  nodes  (  =  2268  degrees  of  freedom),  the  CPU  time 
required  by  AOINA  is  about  the  same  as  that  by  ABAQUS.  Table  6.7  shows 
additional  CPU  detail  of  ADINA  for  analyzing  the  3/0  problem.  Apparently,  the 
requested  workspace  has  some  effect  on  the  CPU  time,  but  not  as  significant  as 
the  2/0  problem.  In  summary,  we  may  draw  the  following  commentary  remarks: 

i)  For  small  size  problems,  ABAQUS  is  not  as  efficient  as  ADINA.  However 

for  large  size  problems,  the  two  codes  seem  to  have  similar  efficiency, 

ii)  The  execution  time  of  ADINA  is  quite  sensitive  to  the  workspace 
allocated.  That  is  to  say,  the  bigger  core  size  is  allocated  to 
ADINA,  the  better  its  computational  efficiency  becomes.  On  the 
contrary,  ABAQUS  is  not  so  sensitive  to  the  size  of  workspace,  so  long 
as  it  is  sufficient  for  the  analysis. 
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Fig.  6.50.  Comparison  of  CPU  Time  Between  ABAQUS  and  ADINA  for  a  3/D  Problem 
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Table  6. 

7  AOINA  Runs  for  the 

3/0  Linear 

Problem 

Problem 

Size 

Workspace 

Number 
of  Blocks 

CPL'  for 

Stiffness 

Formation 

CPU  for 
Triangu- 
larize 

Step-by- 

step 

Total 

Total 

CPU 

2x2x2 

90,000 

1 

13.0sec 

14sec 

1  sec 

47.9  sec 

3x3x3 

90,000 

4 

iminjgsec 

2mi ni7sec 

22  sec 

4min  20.32Sec 

3x3x3 

300,000 

1 

4isec 

imin55sec 

3  sec 

3'^i'i  20.13Sec 

4x4x4 

90,000 

16 

sminsQsec 

13min37sec 

jmin  gsec 

2imin  g^jgsec 

4x4x4 

300,000 

4 

2min48sec 

12'’iini7sec 

imini7SfeC 

17mini4.o8Sec 

5x5x5 

90,000 

66* 

— 

— 

-- 

-- 

5x5x5 

300,000 

11 

45min4sec 

45min8sec 

3min  7sec 

59"'if’30.1isec 

>  *  Run  failed  due  to  too  many  blocks  that  the  structural  stiffness  was 

partitioned. 
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7.  SUMMARY  AND  CONCLUSION 

A  general  purpose  nonlinear  finite  element  program  ABAQUS  was  reviewed 
and  evaluated  with  respect  to  its  documentation,  program  architecture  and 
analysis  capability.  In  addition,  its  operating  characteristics  were  also 
studied  by  running  several  representative  benchmark  problems  and  the  results 
are  given  in  the  advanced  evaluation  section. 

Based  on  our  evaluation  work,  some  specific  points  about  ABAQUS  are 
summarized  as  below: 

1)  ABAQUS  is  supported  by  four  documents:  User's  Manual,  Theory 
Manual,  Systems  Manual  and  Example  Problems  Manual.  Both  the  User's 
Manual  and  Theory  Manual  are  rated  excellent  in  terms  of  their 
clarity  and  information  included.  The  Systems  Manual,  although  it 
has  high  quality,  is  outdated.  The  Example  Problems  Manual  does  not 
have  sufficient  representative  problems  to  demonstrate  different 
aspects  of  the  code. 

2)  Division  of  the  code  into  two  distinct  parts,  the  PRE  and  the  MAIN 
programs,  is  an  efficient  approach  for  data  check  runs  and  quick 
computer  turn  around  time.  Moreover,  the  PRE  program  can  be  run  on  a 
stand-along  mini-  or  micro  -  computer  for  finite  element  model 
preparation  purpose. 

3)  The  coding  structure  is  highly  modularized  so  that  future  extension 


I 


or  modification  becomes  relaively  easy  for  the  developers.  For 
example,  the  material  module  (or  library)  is  completely  independent 
from  the  element  module.  This  is  distinctly  different  from  ADINA. 
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4)  Since  the  code  was  written  by  using  the  data  base  management  and 
separate  workspace  in  the  central  memory  for  processing  data,  all 
data  are  well  protected  from  accidental  over-writing  or  delegation. 

5)  The  Version  3  of  ABAQUS  has  a  fairly  extensive  element  library.  But 
its  material  library  is  somewhat  limited. 

6)  For  loading  (or  time)  increment  control,  the  user  has  two  options  to 
choose  from:  either  direct  stepping  or  automatic  stepping. 

Automatic  stepping  is  extremely  useful  to  start-up  a  new  problem  for 
which  the  user  does  not  have  any  experience  in  its  nonlinear 
behavior.  However,  indiscriminate  use  of  this  option  often  causes 
excessive  CPU  time, 

7)  The  code  is  relatively  machine  independent  (primarily  for  main 
frames).  Any  coding  differences  for  various  machines  were  prepared 
in  a  master  cc'py  of  the  code.  A  special  version  for  a  given  machine 
can  be  obtained  using  text  editor  to  extract  appropriate  coding 
statements. 

8)  The  solution  algorithm  adopted  in  ABAQUS  is  the  full  Newton-Raphson 
method  which  can  handle  a  wide  range  of  nonlinear  problems. 

9)  Logical  input  sequence  with  the  use  of  key-word  cards  and  free  field 
format  is  indeed  one  of  the  desirable  user-friendly  features. 

10)  The  finite  element  formulations  of  ABAQUS  were  consistently  derived 
on  a  firm  theoretical  basis.  All  assumptions  were  clearly  stated  in 
the  theory  manual . 
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11)  By  examing  the  programming  architecture  of  ABAQUS,  the  reviewers 
found  that  it  has  very  high  quality  coding  style  and  programming 
discipline. 

In  additon  to  the  above  comments,  several  areas  may  deserve  some 
attention  for  future  improvement: 

1)  The  nonlinear  solution  algorithm  adopted  in  ABAQUS  is  a  full  Newton- 
Raphson  base.  Other  algorithms,  such  as  modified  and  quasi  Newton- 
Raphson  methods,  should  be  included  in  the  code  for  user's  choice. 

2)  Solution  convergence  during  iterations  is  controlled  by  two  user 
specified  parameters:  PTOL  {force  tolerance)  and  MTOL  (moment 
tolerance).  These  are  the  absolute  physical  quantities  of  force  and 
moment,  respectively.  After  some  extensive  usage  of  the  code,  the 
reviewers  still  feel  uncertain  on  the  question  of  specifying  the 
right  values  of  PTOL  and  MTOL  when  a  new  nonlinear  problem  is 
encountered. 

3)  The  automatic  stepping  method,  although  very  useful  for  starting  up  a 
new  problem,  is  inefficient  when  used  indiscriminately.  A  self- 
adaptive  procedure  with  some  degree  of  intelligence  would  be  most 
desirable. 

4)  More  sophisticated  pre-  and  post-processors  are  needed  for  ABAQUS. 

Or,  as  an  alternative,  its  interface  with  commercial  finite  element 
graphics  packages  would  be  adequate  for  application  purposes. 


7-4 


In  conclusion,  it  is  the  reviewers'  opinion  that  ABAQUS  is  a  strong 
contender  in  the  nonlinear  engineering  analysis  market-place.  The  program  is 
considered  primarily  as  a  production  code,  with  user's  heavy  reliance  on  the 
developers  for  software  support.  Any  coding  changes  or  modifications  by  the 
user  appear  to  be  difficult,  since  the  code  is  not  trannsparent  to  the  common 
users. 
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Appendix  A  -  Dynamic  Response  of  a  Semi-Spherical  Shell 
(Problem  6.7) 


The  solution  for  dynamic  response  of  a  semi-spherical  shell  was 
obtained  by  Kraus  and  Kalnins  [39].  The  non-dimensional i zed  trans¬ 
verse  displacement  is  given  by 


00  w.  (4))  m.F.(T) 

•  -^7- 


(A-1) 


W  =  (D/a^  Pq)  w 


(A-2) 


where 


4)  =  Angular  coordinate,  0  <  4>  <  90" 


T  =  (D/aV)^'^^  •  t 


(A-3) 


12(l-v^) 


(A-4) 


a  =  Radius  of  the  shell 


h  =  Thickness  of  the  shell 


p  =  Density 


t  =  Time 


Pjj  =  Uniform  Pressure 


w  =  Transverse  deflection 


Mode  participation  factor 


A-2 


The  function  F.  in  Eq.  (A-1)  has  the  expression 

F.  =  {T./[K  +  12(1-v2)  (a/h)^]}  * 

[l-(e''''V,.)(A  sin  T.T  +  T.  cos  T.T)] 
and 

=  [K  +  I2(l-v^)  (a/h)^  - 

K  =  (a^/D)  k 

n,.  =  (P/E)^/^  a  0). 

A  =  (a^/phD)^/^  •  X/2 

where 


(A-5) 

(A-6) 

(A-7) 

(A-8) 

(A-9) 


k  =  Elastic  parameter  of  viscoelastic  foundation 
cu^  =  Natural  frequency 

X  =  Viscous  damping  of  the  foundation 

For  the  following  specialized  values  of  parameters: 


h  = 

1  in. 

a 

=  20  in. 

Po  =  1 

E  = 

30  X  10®  Psi 

V 

=  0.3 

2  ‘4 

p  =  1  lb-sec  /m 

n  “ 

0 

k 

=  0 

and  simply  support  condition,  the  transverse  deflection  at  the  corner 
of  the  shell  (Fig.  6.36)  is 


A-3 


c  «>  m . 

w  =  1.33  X  10”^  Z  (1  -  cos  273.86  fi^t)  (A-10) 

i=l  “i  ’ 


The  values  of  m.  and  are  given  as  follows: 


T_ 

Ei 

1 

0.0480 

0.773 

2 

-0.1048 

1.010 

3 

0.3582 

1.242 

4 

-1.5912 

1.578 

5 

2.7978 

1.720 

6 

-0.8794 

2.192 

7 

0.3739 

2.793 

8 

0.2619 

2.941 

9 

-0.4800 

3.733 

10 

0.3596 

4.682 

11 

0.0477 

4.843 

12 

-0.3783 

5.804 

13 

0.0557 

6.876 

14 

0.3063 

7.047 

15 


-0.3411 


8.366 


I 
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Appendix  B  -  Longitudinal  Wave  Propagation  in  a  Bar 
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f 
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iirC, 


The  solution  of  elastic  longitudinal  wave  propagation  can  be  found 

from  Reference  [42].  Consider  a  semi-infinite  bar  subjected  to  an  axial 
force  P(t),  the  displacement  u  is  given  by: 

^  ^  ‘~T~  W~ 

0 

(B-1) 

where 

u  =  Axial  displacement 


^  ^  2L  /  P(t')sin[  2l  (t-f)]  df 


*  Elastic  wave  velocity 

C 


=  /E/p 


A  »  Cross  sectional  area 


L  *  Length  of  the  bar 
P  *  Pressure  pulse  applied  at  one  end 


For  Problem  6.8,  we  consider  the  following  pressure  history 
p  =  Pq  •  t/T^  ,  te  [o,  Tj] 

where 

Pjj  »  Peak  pressure  at  t  = 


B-2 


The  stress  pulse  has  the  expression 


-  F,(t)cos^ 


(B-3) 


where 


F,  «  - 


sin  p.t 


t  e  [o.T^] 


T, 


sin  P^.  t 

Pi  h 


^2  (t-Tj) 


T2  sin  (t-T^ 
'’i’l 


t  e  [tj,  T2] 


Sin  P^t  ^2 

"W  '  “ 


T,  Sin  P^  (t-Tj) 


Sin  P.  (t-T2) 

npTTfpfTT” 


,  t  ^  T2 


(B-4) 


The  above  series  converges  very  slowly;  it  requires  about  200  natural 
modes. 


B-3 


In  the  case  of  plastic  wave  propagation  (for  elastic-perfectly  plastic 
material),  let  the  stress  pulse  be  denoted  by  (also  shown  in  Fig.  B-1). 


opront  = 
^Back  = 


(B-5) 


where  t  is  the  time  at  which  the  stress  pulse  is  applied;  kp  and  k^  are  the 
slopes  of  the  front  and  back  of  the  pulse.  The  stress  profile  travelling  in 
the  bar  at  different  instants  varies  due  to  the  difference  in  elastic  wave 
speed  Cp  and  plastic  wave  speed  Cp.  To  determine  the  stress  profiles,  we 
consider  the  intersection  of  the  plastic  front  and  elastic  back  (corresponding 
to  plastic  unloading).  The  intersection  has  a  stress  value  'a  and  is  located 
at  a  distance  "d  from  the  free  end  of  the  bar.  Both  7  and  d  can  be  expressed 
by 

0  »  kp  xp  +  ap  =  kg  Tg  +  og 

(B-6) 

■J  »  (t  -  xp)  Cp  =  (t  -  Tg)  Cp 


where  xp  and  xg  are  the  starting  time  of  the  stress  pulse  a  on  the  front  and 
back,  respectively  as  shown  in  Fig.  B-1.  From  Eqs.  (B-5)  and  (B-6),  solving 
for  xp  and  xg 


(1  - 


1-  IT  i 


'^F 


(B-7) 


(b)  Intersection  of  the  Front  and  Back  of  Stress  Pulse 


M 
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Appendix  C  -  Additional  Analysis  Capability  of  Version 
4  -  ABAQUS.* 

For  reader's  information,  additional  analysis  capabilities  of 
version  4.5  ABAQUS  over  Version  3  are  listed  below: 

Elements 

-  Continuum  elements  with  finite  strains 

-  Incompressible  hybrid  elements  for  plane  strain,  axi symmetric 
and  three-dimensional  analysis. 

-  Axi symmetric  shells  with  finite  strains 

-  Line  Spring  elements  for  modeling  part-through  cracks  in  shells 

-  Shell  elements  for  heat  —  transfer  analysis 

-  2/D  and  3/D  elements  with  pore  pressures  for  consideration 
analysi s 

-  Interface  elements  (2/D  and  3/D)  for  contact/friction  stress 
problems 

-  Interface  elements  (2/D  and  3/D)  for  heat  transfer  analysis 

-  Elements  for  coupled  thermal -stress  analysis  (truss,  2/D,  3/D 
and  shell  elements) 

-  Pipe  elements  with  the  effect  of  internal  pressure 
Materials 

-  A  hypoelastic  model  for  soils 

-  Modified  Cam  Clay  model 

-  ORNL  creep  and  plasticity  model  for  types  306  and  316  stain¬ 
less  steel 


C-2 


-  Coulomb  friction  model  for  interface  elements 

-  No-tension  or  no-compression  material  model 

-  Gap  conductance  and  radiation  for  interface  elements 

-  Latent  heat  model 

-  Permeability  model  (isotropic,  orthotropic  or  anistropic) 
for  soils 

-  Chen  and  Chen  concrete  plasticity  model 

-  hyperelastic  model  for  elastomer 

Procedures 

-  A  modified  Riks  algorithm  for  structural  instability  analysis 

-  Consolidation  analysis  for  soils 

-  Coupled  temperature/displacement  analysis 

-  Fluid-solid  interaction  analysis 

-  Superelement/substructuring  capability 

-  J  integral  calculations  with  Park's  node  moving  technique 


*  The  above  information  was  supplied  by  Hibbitt,  Karlsson  &  Sorensen,  Inc. 


(March  30,  1983). 
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